Changeset 39204

Timestamp:

Dec 1, 2015, 7:15:51 AM (11 years ago)

Author:

eugene

Message:

add monte-carlo distance

File:

• branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh (modified) (4 diffs)

branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh

@@ -35 +35 @@
 # generate the predicted proper motions
 
+macro mkstars_dmerr
+  if ($0 != 5)
+    echo "USAGE: mkstars (Nstars) (Dm) (Dm_err) (Nsample)"
+    break
+  end
+
+  create seq 0 $1
+  set phi   = 2*rnd(seq) - 1
+
+  set d_kpc = ten(0.2*$2 - 2)
+
+  set L = 360*rnd(seq)
+  set B = dasin(phi)
+
+  set R = L
+  set D = B
+  csystem G C R D
+
+  set uL_gal =     ($A_oort * dcos(2.0*L) + $B_oort) * dcos(B*1.0) * $iFkap;
+  set uB_gal = -0.5*$A_oort * dsin(2.0*L) *            dsin(B*2.0) * $iFkap;
+
+  set uL_sol =  ($U_sol * dsin(L) - $V_sol * dcos(L))                           * $iFkap / $d_kpc;
+  set uB_sol = (($U_sol * dcos(L) + $V_sol * dsin(L))*dsin(B) - $W_sol*dcos(B)) * $iFkap / $d_kpc;
+
+  set uL_o = uL_gal + uL_sol
+  set uB_o = uB_gal + uB_sol
+
+  set uL_sol_o = uL_sol
+  set uB_sol_o = uB_sol
+
+  # create a buffer to store the results : results for each star are in the y-dir, 
+  mcreate resultsL $1 $4
+  mcreate resultsB $1 $4
+
+  # now MC a number of distances
+  for i 0 $4
+    echo $i
+    gaussdev Dm_x $1 $2 $3
+    set d_kpc = ten(0.2*Dm_x - 2)
+    
+    set uL_sol =  ($U_sol * dsin(L) - $V_sol * dcos(L))                           * $iFkap / d_kpc;
+    set uB_sol = (($U_sol * dcos(L) + $V_sol * dsin(L))*dsin(B) - $W_sol*dcos(B)) * $iFkap / d_kpc;
+
+    set duL = uL_sol - uL_sol_o
+    set duB = uB_sol - uB_sol_o
+
+    mset resultsL duL -x $i    
+    mset resultsB duB -x $i    
+  end
+    
+  set duL = zero(uL_o)
+  set duB = zero(uB_o)
+  set uLoff = zero(uL_o)
+  set uBoff = zero(uB_o)
+
+  for i 0 $1
+    mget resultsL tmp -y $i
+    vstat -q tmp
+    
+    uLoff[$i] = $MEAN
+    duL[$i] = $SIGMA
+
+    mget resultsB tmp -y $i
+    vstat -q tmp
+    
+    uBoff[$i] = $MEAN
+    duB[$i] = $SIGMA  
+  end
+end
+
 macro mkstars
-  if ($0 != 4)
-    echo "USAGE: mkstars (Nstars) (d kpc) (scale)"
+  if ($0 != 3)
+    echo "USAGE: mkstars (Nstars) (d kpc)"
     break
   end
@@ -46 +116 @@
   set L = 360*rnd(seq)
   set B = dasin(phi)
-
-  resize -n galactic 1600 800
-  region 0 0 90 ait; cgrid -c red
-
-  resize -n celestial 1600 800
-  region 0 0 90 ait; cgrid -c red
 
   set R = L
@@ -66 +130 @@
   set uB = uB_gal + uB_sol
 
-  dev -n galactic
-  cneedles L B uL uB -scale $3 -c black -lw 0
-  cplot L B -pt 0 -sz 0.3 -c black
-  label +x Galactic -fn courier 24
-
   # C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf
   set C1 =  dcos(D)*$trans_cos_phi + dsin(D)*(dcos(R)*$trans_sin_phi_sin_Xo - dsin(R)*$trans_sin_phi_cos_Xo);
@@ -79 +138 @@
   set uR = cosBinv * (C1 * uL - C2 * uB);
   set uD = cosBinv * (C1 * uB + C2 * uL);
+end
 
-  dev -n celestial
+macro pltstars
+  if ($0 != 2)
+    echo "USAGE: pltstars (scale)"
+    break
+  end
+
+  resize -n galactic 1600 800
+  region 0 0 90 ait; cgrid -c red
+
+  cneedles L B uL uB -scale $1 -c black -lw 0
+  cplot L B -pt 0 -sz 0.3 -c black
+  label +x Galactic -fn courier 24
+
+  resize -n celestial 1600 800
+  region 0 0 90 ait; cgrid -c red
+
   style -c red -lw 0 -pt 0 -sz 0.3; cggrid; style -c blue; galactic
   cplot R D -pt 0 -sz 0.3 -c black
-  cneedles R D uR uD -scale $3 -c black -lw 0
+  cneedles R D uR uD -scale $1 -c black -lw 0
   label +x Celestial -fn courier 24
-
-  
 end