Index: /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh	(revision 39207)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh	(revision 39208)
@@ -4,5 +4,5 @@
 
 ## FEAST-HIP
-if (1)
+if (0)
   $A_oort = +14.82; # km/sec/kpc
   $B_oort = -12.37; # km/sec/kpc
@@ -13,5 +13,5 @@
 
 ## ROESSER:
-if (0)
+if (1)
   $A_oort = +14.50; # km/sec/kpc
   $B_oort = -13.00; # km/sec/kpc
@@ -37,6 +37,64 @@
 $trans_sin_Xo 	      = dsin($trans_Xo);
 
+$disk_radius = 10; # 10 kpc radius disk
+$disk_scale = 0.5; # 200 pc scale height
+
+$halo_scale = 2.0; # 2000 pc scale height
+
 # generate a collection of stars in R,D [L,B] with a range of distances
 # generate the predicted proper motions
+
+macro mkstars_halo
+  if ($0 != 4)
+    echo "USAGE: mkstars_disk (Nstars) (Dm_err) (Nsample)"
+    break
+  end
+
+  local Nstars dDm_o Nsample
+
+  $Nstars = $1
+  $dDm_o = $2
+  $Nsample = $3
+
+  # halo stars are uniformly sampled in a gaussian halo
+  gaussdev Dhalo $Nstars 0.0 $halo_scale
+  set Dm = 5*log(Dhalo * 100)
+  
+  create seq 0 $Nstars
+  set phi   = 2*rnd(seq) - 1
+
+  set L = 360*rnd(seq)
+  set B = dasin(phi)
+end
+
+macro mkstars_disk
+  if ($0 != 5)
+    echo "USAGE: mkstars_disk (Nstars) (Dm_err) (Nsample) (minDm)"
+    break
+  end
+
+  local Nstars dDm_o Nsample minDm
+
+  $Nstars = $1
+  $dDm_o = $2
+  $Nsample = $3
+  $minDm = $4
+
+  gaussdev Zgal $Nstars 0.0 $disk_scale
+  set rgal = $disk_radius * rnd(Zgal)
+
+  set B = atan2(Zgal , rgal)
+  set L = 360*rnd(Zgal)
+
+  set d_kpc_o = sqrt(rgal^2 + Zgal^2)
+
+  set Dm = 5*log(d_kpc_o * 100)
+
+  subset Bs  = B  if (Dm > $minDm)
+  subset Ls  = L  if (Dm > $minDm)
+  subset Dms = Dm if (Dm > $minDm)
+
+  mkmotion_dmerr Ls Bs Dms $dDm_o $Nsample
+end
 
 macro mkstars_uniform
@@ -61,13 +119,53 @@
   set B = dasin(phi)
 
-  set R = L
-  set D = B
+  mkmotion_dmerr L B Dm $3 $4
+end
+
+macro mkstars_single_dm
+  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) + zero(seq)
+  set Dm = 5*log(d_kpc * 100)
+
+  set L = 360*rnd(seq)
+  set B = dasin(phi)
+
+  mkmotion_dmerr L B Dm $3 $4
+end
+
+# given L,B,Dm and a Dm_err value, find duL, duB
+macro mkmotion_dmerr
+  if ($0 != 6)
+    echo "USAGE: mkmotion (L) (B) (Dm) (Dm_err) (Nsample)"
+    break
+  end
+
+  local i Dm_err Nstars 
+
+  set  _L = $1
+  set  _B = $2
+  set _Dm = $3
+  $Dm_err = $4
+  $Nsample = $5
+
+  $Nstars = _L[]
+
+  set d_kpc_o = ten(0.2*_Dm - 2)
+
+  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_o;
-  set uB_sol = (($U_sol * dcos(L) + $V_sol * dsin(L))*dsin(B) - $W_sol*dcos(B)) * $iFkap / d_kpc_o;
+  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_o;
+  set uB_sol = (($U_sol * dcos(_L) + $V_sol * dsin(_L))*dsin(_B) - $W_sol*dcos(_B)) * $iFkap / d_kpc_o;
 
   set uL_o = uL_gal + uL_sol
@@ -78,18 +176,15 @@
 
   # create a buffer to store the results : results for each star are in the y-dir, 
-  mcreate resultsL $Nstars $4
-  mcreate resultsB $Nstars $4
+  mcreate resultsL $Nstars $Nsample
+  mcreate resultsB $Nstars $Nsample
 
   # now MC a number of distances
-  for i 0 $4
-    echo $i
-
-    gaussdev dDm $Nstars 0.0 $dDm_o
-    set Dm_x = Dm + dDm
-
+  for i 0 $Nsample
+    gaussdev dDm $Nstars 0.0 $Dm_err
+    set Dm_x = _Dm + dDm
     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 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
@@ -120,74 +215,4 @@
 end
 
-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 != 3)