Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c	(revision 37460)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c	(revision 37460)
@@ -0,0 +1,83 @@
+# include "fakeastro.h"
+
+// things to figure out:
+// * which filter am I (image.photcode)
+// * lookup table of sky counts / pixel or seeing disk
+// * what about QSOs?
+
+Stars *make_fake_stars_catalog (Catalog *catalog, Image *image, int *nfakeStars) {
+
+  float Mtime = 2.5*log10(image->exptime);
+  float ZP = image->Mcal + zeropt;
+
+  // generate a set of measurements for each star entry
+
+  Average *average = catalog->average;
+  StarPar *starpar = catalog->starpar;
+
+  for (i = 0; i < catalog->Naverage; i++) {
+
+    int nStar = average[i].starparOffset
+
+    InitStar (&stars[i]);
+
+    // which filter?
+    double Minst = secfilt[i*Nsecfilt + Ns].M  - Mtime - ZP;
+    double Counts = pow(10.0, -0.4*Minst);
+    double SkyCts = Something;
+
+    double SN = Counts / sqrt(SkyCts + Counts);
+
+    // true position from src catalog
+    double Rtru = average[i].R;
+    double Dtru = average[i].D;
+
+    // observed position is scattered from true position by:
+    // * proper motion
+    // * gaussian scatter (~ seeing) 
+    double uR = starpar[nStar].uR;
+    double uD = starpar[nStar].uD;
+    
+    double Toffset = average - image.tzero;
+
+    // uR,uD in linear (arcsec / yr)
+    double dRoff = uR*Toffset;
+    double dDoff = uD*Toffset;
+
+    // uR,uD in linear arcsec
+    double dRsee = rnd_gauss (0.0, 1.0 / SN);
+    double dDsee = rnd_gauss (0.0, 1.0 / SN);
+
+    double Robs = Rtru + dRoff / cos(Dtru*DEG_RAD);
+    double Dobs = Dtru + dDoff;
+
+    double X, Y;
+    RD_to_XY (&X, &Y, Robs, Dobs, image->coords);
+
+    stars[i].measure.Xccd       = X;
+    stars[i].measure.Yccd       = Y;
+    stars[i].measure.dXccd      = 1.0 / SN / plateScale;
+    stars[i].measure.dYccd      = 1.0 / SN / plateScale;
+
+    // stars[i].measure.posangle   = ToShortDegrees(ps1data[i].posangle);
+    // stars[i].measure.pltscale   = ps1data[i].pltscale;
+
+    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
+      stars[i].measure.M      = NAN;
+    } else {
+      stars[i].measure.M      = Minst + ZeroPt;
+    }
+    stars[i].measure.dM         = 1.0 / SN;
+
+    // stars[i].measure.dMcal      = ps1data[i].dMcal;
+    stars[i].measure.Sky        = X?;
+    stars[i].measure.dSky       = X?;
+                        
+    stars[i].measure.photFlags  = ps1data[i].flags;
+    stars[i].measure.photFlags2 = ps1data[i].flags2;
+
+    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
+    stars[i].measure.detID      = ps1data[i].detID; 
+  }
+
+}