Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/Makefile
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/Makefile	(revision 39735)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/Makefile	(revision 39736)
@@ -53,5 +53,7 @@
 $(SRC)/match_fake_stars.$(ARCH).o \
 $(SRC)/fakeastro_2mass.$(ARCH).o \
+$(SRC)/fakeastro_gaia.$(ARCH).o \
 $(SRC)/make_2mass_measures.$(ARCH).o \
+$(SRC)/make_gaia_measures.$(ARCH).o \
 $(SRC)/remote_hosts.$(ARCH).o 
 
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/include/fakeastro.h
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/include/fakeastro.h	(revision 39735)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/include/fakeastro.h	(revision 39736)
@@ -8,5 +8,5 @@
 # define RESETTIME { gettimeofday (&startTimer, (void *) NULL); }
 
-typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES, OP_2MASS} FakeastroOp;
+typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES, OP_2MASS, OP_GAIA} FakeastroOp;
 
 typedef struct {
@@ -101,7 +101,9 @@
 char   FAKEASTRO_REF_EPOCH[80];
 char   FAKEASTRO_2MASS_EPOCH[80];
+char   FAKEASTRO_GAIA_EPOCH[80];
 
 float  RADIUS;
 float  MAX_MAG_2MASS;
+float  MAX_MAG_GAIA;
 
 SkyRegion UserPatch;
@@ -198,2 +200,5 @@
 int fakeastro_2mass ();
 int make_2mass_measures (Catalog *catalog);
+
+int fakeastro_gaia ();
+int make_gaia_measures (Catalog *catalog);
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/src/ConfigInit.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/src/ConfigInit.c	(revision 39735)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/src/ConfigInit.c	(revision 39736)
@@ -19,5 +19,5 @@
   // if (!ScanConfig (config, "ADDSTAR_RADIUS",         "%lf", 0, &ADDSTAR_RADIUS))   ADDSTAR_RADIUS = 1.0;
 
-  if ((FAKEASTRO_OP == OP_GALAXY) || (FAKEASTRO_OP == OP_2MASS)) {
+  if ((FAKEASTRO_OP == OP_GALAXY) || (FAKEASTRO_OP == OP_2MASS) || (FAKEASTRO_OP == OP_GAIA)) {
     // force CATDIR to be absolute (so parallel mode will work)
     GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
@@ -60,4 +60,7 @@
     strcpy (FAKEASTRO_2MASS_EPOCH, "2000/01/01,00:00:00"); // epoch of 2MASS astrometry
   }
+  if (!ScanConfig (config, "FAKEASTRO_GAIA_EPOCH", "%s", 0, FAKEASTRO_GAIA_EPOCH)) {
+    strcpy (FAKEASTRO_GAIA_EPOCH, "2015/01/01,00:00:00"); // epoch of GAIA astrometry
+  }
 
   /* set the default search radius */
@@ -76,5 +79,5 @@
 
   // OP_2MASS is adding detections to an existing db, the others require and empty db
-  if (FAKEASTRO_OP != OP_2MASS) {
+  if ((FAKEASTRO_OP != OP_2MASS) && (FAKEASTRO_OP != OP_GAIA)) {
     // check for existence of CATDIR
     struct stat filestat;
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/src/args.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/src/args.c	(revision 39735)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/src/args.c	(revision 39736)
@@ -24,4 +24,9 @@
     remove_argument (N, argc, argv);
     FAKEASTRO_OP = OP_2MASS;
+  }
+
+  if ((N = get_argument (*argc, argv, "-gaia"))) {
+    remove_argument (N, argc, argv);
+    FAKEASTRO_OP = OP_GAIA;
   }
 
@@ -114,4 +119,11 @@
   }
 
+  MAX_MAG_GAIA = 16.0;
+  if ((N = get_argument (*argc, argv, "-gaia-limit"))) {
+    remove_argument (N, argc, argv);
+    MAX_MAG_GAIA = atof(argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
   FORCE = FALSE;
   if ((N = get_argument (*argc, argv, "-force"))) {
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro.c	(revision 39735)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro.c	(revision 39736)
@@ -41,4 +41,10 @@
       exit (0);
 
+    case OP_GAIA:
+      fakeastro_gaia ();
+      /* make_gaia_measures()
+       */
+      exit (0);
+
     default:
       fprintf (stderr, "impossible!\n");
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro_gaia.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro_gaia.c	(revision 39736)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/src/fakeastro_gaia.c	(revision 39736)
@@ -0,0 +1,57 @@
+# include "fakeastro.h"
+
+int fakeastro_gaia () {
+
+  INITTIME;
+
+  SkyTable *skyTable = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (skyTable, CATDIR, "cpt");
+
+  SkyList *skylist  = SkyListByPatch (skyTable, -1, &UserPatch);
+
+  Catalog catalog;
+
+  // load stars from database in these regions
+  int i;
+  for (i = 0; i < skylist->Nregions; i++) {
+    dvo_catalog_init (&catalog, TRUE);
+
+    // set the parameters which guide catalog open/load/create
+    dvo_catalog_init (&catalog, TRUE);
+    catalog.filename  = skylist[0].filename[i];
+    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
+    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
+    catalog.catflags  = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT | DVO_LOAD_MEASURE | DVO_LOAD_STARPAR;
+    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
+      exit (1);
+    }
+    if (!catalog.Naverage_disk) {
+	if (VERBOSE2) { fprintf (stderr, "no data in %s, skipping\n", catalog.filename); }
+	dvo_catalog_unlock (&catalog);
+	dvo_catalog_free (&catalog);
+	continue;
+    }
+
+    make_gaia_measures (&catalog);
+
+    SetProtect (TRUE);
+    if (!dvo_catalog_update (&catalog, VERBOSE)) {
+      fprintf (stderr, "ERROR: failure to update %s\n", catalog.filename);
+      exit (3);
+    }
+    SetProtect (FALSE);
+
+    if (!dvo_catalog_unlock (&catalog)) {
+      fprintf (stderr, "ERROR: failure to unlock %s\n", catalog.filename);
+      exit (2);
+    }
+    dvo_catalog_free (&catalog);
+    if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; 
+  }
+
+  exit (0);
+}
+
+
Index: /branches/czw_branch/20160809/Ohana/src/fakeastro/src/make_gaia_measures.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/fakeastro/src/make_gaia_measures.c	(revision 39736)
+++ /branches/czw_branch/20160809/Ohana/src/fakeastro/src/make_gaia_measures.c	(revision 39736)
@@ -0,0 +1,123 @@
+# include "fakeastro.h"
+
+// region corresponds to the catalog
+int make_gaia_measures (Catalog *catalog) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  Average *average = catalog->average;
+  SecFilt *secfilt = catalog->secfilt;
+  StarPar *starpar = catalog->starpar;
+  Measure *measure = catalog->measure;
+  
+  // we have starpar.R,D, which represent the true positions at the reference epoch, FAKEASTRO_REF_EPOCH
+  // the GAIA stars are generated at locations for the GAIA EPOCH
+  time_t timeRef = ohana_date_to_sec(FAKEASTRO_REF_EPOCH);
+  time_t tzero_gaia = ohana_date_to_sec(FAKEASTRO_GAIA_EPOCH);
+
+  // tzero, timeRef are in UNIX seconds, Toffset should be in years
+  double Toffset = (tzero_gaia - timeRef) / 365.25 / 86400.0;
+
+  // use photcode to get zero point
+  PhotCode *codeG = GetPhotcodebyName ("GAIA_G_DR1");
+
+  float ZP = MAX_MAG_GAIA + 6.5; // this is chosen to give stars at the limit a flux of 400.0 counts (crude, yes)
+  float SkyCts = 1000.0; // this is chosen to make SN @ limit ~ 10.0 (
+
+  int Nmeasure = catalog->Nmeasure;
+  int NMEASURE = catalog->Nmeasure;
+
+  int Nsec = 2; // i-band
+
+  int i;
+  for (i = 0; i < catalog->Naverage; i++) {
+
+    if (Nmeasure >= NMEASURE) {
+      NMEASURE = Nmeasure + 1000;
+      REALLOCATE (catalog[0].measure, Measure, NMEASURE);
+      measure = catalog[0].measure;
+    }
+
+    int Nstarpar = average[i].starparOffset;
+
+    // make a crude G-i color for now:
+    double G_PS1 = secfilt[i*Nsecfilt + Nsec].M + 1.0; // make all stars have G-i = 1.0
+    if (G_PS1 > MAX_MAG_GAIA) continue; // only generate GAIA detections for objects with G_PS1 < 21.0
+
+    if (isnan(secfilt[i*Nsecfilt + Nsec].M)) {
+      // look for a non-NAN secfilt mag and just use that (it is not super important)
+      int ns;
+      for (ns = 0; ns < Nsecfilt; ns++) {
+	if (!isnan(secfilt[i*Nsecfilt + ns].M)) break;
+      }
+      if (ns == Nsecfilt) continue; // no non-nan
+
+      G_PS1 = secfilt[i*Nsecfilt + ns].M; // pretend secfilt.M = G
+      if (G_PS1 > MAX_MAG_GAIA) continue; // only generate GAIA detections for objects with G_PS1 < 21.0
+    }
+
+    double Minst = G_PS1 - ZP; 
+
+    double Counts = pow(10.0, -0.4*Minst);
+
+    double SN = Counts / sqrt(SkyCts + Counts);
+
+    // true position from src catalog
+    double Rtru = starpar[Nstarpar].R;
+    double Dtru = starpar[Nstarpar].D;
+
+    // observed position is scattered from true position by:
+    // * proper motion
+    // * gaussian scatter (~ seeing) 
+    double uR = starpar[Nstarpar].uRA; // starpar are in arcsec / year
+    double uD = starpar[Nstarpar].uDEC;
+    
+    // uR,uD in linear (arcsec / yr)
+    double dRpm = uR*Toffset;
+    double dDpm = uD*Toffset;
+
+    // uR,uD in linear arcsec
+    double dRsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
+    double dDsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
+
+    double dRoff = (dRpm + dRsee) / 3600.0;
+    double dDoff = (dDpm + dDsee) / 3600.0;
+
+    double Robs = Rtru + dRoff / cos(Dtru*RAD_DEG);
+    double Dobs = Dtru + dDoff;
+
+    dvo_measure_init (&measure[Nmeasure]);
+
+    measure[Nmeasure].R = Robs;
+    measure[Nmeasure].D = Dobs;
+
+    measure[Nmeasure].M      = G_PS1;
+    measure[Nmeasure].dM     = 1.0 / SN;
+
+    measure[Nmeasure].Sky        = SkyCts;
+    measure[Nmeasure].dSky       = sqrt(SkyCts);
+    measure[Nmeasure].photFlags  = 0;
+    measure[Nmeasure].photFlags2 = 0;
+    measure[Nmeasure].airmass = 1.0;
+    measure[Nmeasure].az      = 0.0; // irrelevant
+    measure[Nmeasure].Mcal    = 0.0;
+    measure[Nmeasure].t       = tzero_gaia;
+    measure[Nmeasure].dt      = 0.0;
+    measure[Nmeasure].photcode = codeG->code;
+
+    measure[Nmeasure].averef   = i;
+    measure[Nmeasure].objID    = average[i].objID;
+    measure[Nmeasure].catID    = average[i].catID;
+
+    measure[Nmeasure].imageID = 0;
+
+    // This is may optionally be replaced by the internal sequence (see FilterStars.c)
+    measure[Nmeasure].detID      = 0;
+    Nmeasure ++;
+  }
+
+  if (VERBOSE) fprintf (stderr, "added %d gaia entries to %d stars\n",  (int) Nmeasure, (int) catalog->Naverage);
+
+  catalog->Nmeasure = Nmeasure;
+  return TRUE;
+}
