Index: /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh	(revision 39195)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/uniphot/doc/galmotion.sh	(revision 39195)
@@ -0,0 +1,77 @@
+
+# global constants:
+$iFkap    = 0.001 / 4.74047; # (km/sec/kpc) / (arcsec/ year)
+
+## FEAST-HIP
+$A_oort = +14.82; # km/sec/kpc
+$B_oort = -12.37; # km/sec/kpc
+$U_sol  =   9.32; # km/sec
+$V_sol  =  11.18; # km/sec
+$W_sol  =   7.61; # km/sec
+
+## ROESSER:
+$A_oort = +14.50; # km/sec/kpc
+$B_oort = -13.00; # km/sec/kpc
+$U_sol  =   9.44; # km/sec
+$V_sol  =  11.90; # km/sec
+$W_sol  =   7.20; # km/sec
+
+$trans_phi = -62.8717488056*RAD_DEG;
+$trans_Xo  = 282.8594812080*RAD_DEG;
+$trans_xo  =  32.9319185700;
+
+$trans_sin_phi_cos_Xo = dsin($trans_phi)*dcos($trans_Xo);
+$trans_sin_phi_sin_Xo = dsin($trans_phi)*dsin($trans_Xo);
+$trans_cos_phi        = dcos($trans_phi);
+
+$trans_cos_phi_cos_Xo = dcos($trans_phi)*dcos($trans_Xo);
+$trans_cos_phi_sin_Xo = dcos($trans_phi)*dsin($trans_Xo);
+$trans_sin_phi        = dsin($trans_phi);
+
+$trans_cos_Xo 	      = dcos($trans_Xo);
+$trans_sin_Xo 	      = dsin($trans_Xo);
+
+# generate a collection of stars in R,D [L,B] with a range of distances
+# generate the predicted proper motions
+
+macro mkstars
+  if ($0 != 4)
+    echo "USAGE: mkstars (Nstars) (d kpc) (scale)"
+    break
+  end
+
+  create seq 0 $1
+  set phi   = 2*rnd(seq) - 1
+
+  set L = 360*rnd(seq)
+  set B = dasin(phi)
+
+  region 0 0 90 ait; cgrid -c red
+  # cplot L B
+  # break
+
+  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 / $2;
+  set uB_sol = (($U_sol * dcos(L) + $V_sol * dsin(L))*dsin(B) - $W_sol*dcos(B)) * $iFkap / $2;
+
+  set uL = uL_gal + uL_sol
+  set uB = uB_gal + uB_sol
+
+  cneedles L B uL uB -scale $3
+  break
+
+  # 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);
+  set C2 =                                -1*(dcos(R)*$trans_sin_phi_cos_Xo + dsin(R)*$trans_sin_phi_sin_Xo);
+
+  set cosBinv = 1.0 / sqrt(C1*C1 + C2*C2);
+
+  set uR = cosBinv * (C1 * uL - C2 * uB);
+  set uD = cosBinv * (C1 * uB + C2 * uL);
+end