# include <dvo.h>

# define NCTERMS 4
# define F_PS 0.001
# define NO_MAG_PS 100.0

static PhotCodeData *photcodes = NULL;
static double ZERO_POINT;
static short int Nseclist[0x10000];

/* static short int iZERO_POINT; */

void SetZeroPoint (double ZP) {
  ZERO_POINT = ZP;
  /* iZERO_POINT = 1000 * ZP; */
}

int LoadPhotcodes (char *filename) {
  
  FILE *f;
  int i, Ns, Np, NPHOTCODE, Npri, Nfield;
  PhotCode *photcode;
  int code;
  char *c;
  char line[256];
  char name[32], type[32], Zero[32], Airmass[32], Offset[32],
    C1[32], C2[32], Slope[32], Color[32], Primary[32];
  int c1, c2;

  /* allocate space to photcode table, free existing data */
  if (photcodes == NULL) {
    ALLOCATE (photcodes, PhotCodeData, 1);
    photcodes[0].code = NULL;
  }
  if (photcodes[0].code != NULL) free (photcodes[0].code);

  f = fopen (filename, "r");
  if (f == (FILE *) NULL) {
    photcodes[0].Ncode    = 0;
    photcodes[0].Nsecfilt = 0;
    photcodes[0].code     = (PhotCode *) NULL;
    return (FALSE);
  }

  Np = 0;
  NPHOTCODE = 10;
  ALLOCATE (photcode, PhotCode, NPHOTCODE);
  photcodes[0].Nsecfilt = 0;
  for (i = 0; i < 0x10000; i++) {
    photcodes[0].hashcode[i] = -1;
    photcodes[0].hashNsec[i] = -1;
  }

  while (scan_line (f, line) != EOF) {
    for (c = line; isspace (*c); c++);
    if (*c == '#') continue;
    Nfield = sscanf (c, "%d %s %s %s %s %s %s %s %s %s %s", 
		     &code, name, type, Zero, Airmass, Offset, C1, C2, Slope, Color, Primary);
    if (Nfield != 11) { continue; }
    
    c1 = atof (C1);
    c2 = atof (C2);
    if (!strcmp (C1, "-")) { c1 = 0; }
    if (!strcmp (C2, "-")) { c2 = 0; }

    photcode[Np].type = 0;
    photcode[Np].code = code;
    strcpy (photcode[Np].name, name);
    if (!strncasecmp (type, "pri", 3)) {
      photcode[Np].type  = PHOT_PRI;
      photcode[Np].C     = 1000*atof (Zero);
      photcode[Np].K     = atof (Airmass);
      photcode[Np].dC    = 1000*atof (Offset);
      photcode[Np].dX    = 1000*atof (Color);
      photcode[Np].c1    = c1;
      photcode[Np].c2    = c2;
      photcode[Np].equiv = atoi (Primary);
      ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
      Nseclist[0] = Np;
    }      
    if (!strncasecmp (type, "sec", 3)) {
      photcode[Np].type  = PHOT_SEC;
      photcode[Np].C     = 1000*atof (Zero);
      photcode[Np].K     = atof (Airmass);
      photcode[Np].dC    = 1000*atof (Offset);
      photcode[Np].dX    = 1000*atof (Color);
      photcode[Np].c1    = c1;
      photcode[Np].c2    = c2;
      photcode[Np].equiv = atoi (Primary);
      photcodes[0].Nsecfilt ++;
      ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
      Nseclist[photcodes[0].Nsecfilt] = Np;
    }      
    if (!strncasecmp (type, "dep", 3)) {
      photcode[Np].type  = PHOT_DEP;
      photcode[Np].C     = 1000*atof (Zero);    /* zero point in millimags */
      photcode[Np].K     = atof (Airmass);      /* airmass coeff (millimag / millimag) */
      photcode[Np].dC    = 1000*atof (Offset);  /* color ref z.p. (millimag) */
      photcode[Np].dX    = 1000*atof (Color);   /* average color (millimag) */
      photcode[Np].c1    = c1;
      photcode[Np].c2    = c2;
      photcode[Np].equiv = atoi (Primary);
      ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    }      
    if (!strncasecmp (type, "ref", 3)) {
      photcode[Np].type  = PHOT_REF;
      photcode[Np].C     = 0;
      photcode[Np].K     = 0;
      photcode[Np].dC    = 0;
      photcode[Np].dX    = 0;
      photcode[Np].c1    = 0;
      photcode[Np].c2    = 0;
      photcode[Np].equiv = atoi (Primary);
      photcode[Np].X[0]  = 0;
      photcode[Np].Nc    = 0;
    }

    /* alt photcodes are a little different: they have the SAME photcode as an existing
       pri/sec photcode, but define an alternate transformation for that code */
    if (!strncasecmp (type, "alt", 3)) {
      photcode[Np].type  = PHOT_ALT;
      photcode[Np].C     = 1000*atof (Zero);    /* zero point in millimags */
      photcode[Np].K     = atof (Airmass);      /* airmass coeff (millimag / millimag) */
      photcode[Np].dC    = 1000*atof (Offset);  /* color ref z.p. (millimag) */
      photcode[Np].dX    = 1000*atof (Color);   /* average color (millimag) */
      photcode[Np].c1    = c1;
      photcode[Np].c2    = c2;
      photcode[Np].equiv = atoi (Primary);
      ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    }
    if (!photcode[Np].type) {
      fprintf (stderr, "error in Photfile: unknown type %s\n", type);
    }

    Np++;
    if (Np == NPHOTCODE) {
      NPHOTCODE += 10;
      REALLOCATE (photcode, PhotCode, NPHOTCODE);
    }
  }
  fclose (f);
  
  /* set up hashcode for photcode refs:
   * the hashcode gives the structure sequence for a given photcode:
   * photcode[i].hashcode[photcode[i].code] == i
   */

  Ns = 0;
  for (i = 0; i < Np; i++) {
    if (photcode[i].type == PHOT_ALT) continue; /* no hashcode for ALT codes */
    if (photcodes[0].hashcode[photcode[i].code] != -1) {
      fprintf (stderr, "duplicate photcodes in file\n");
      code = photcodes[0].hashcode[photcode[i].code];
      fprintf (stderr, "conflict between %s (%d) and %s (%d)\n",
	       photcode[i].name, photcode[i].code, photcode[code].name, photcode[code].code);
      free (photcode);
      return (FALSE);
    }
    photcodes[0].hashcode[photcode[i].code] = i;
    if (photcode[i].type == PHOT_SEC) {
      photcodes[0].hashNsec[photcode[i].code] = Ns;
      Ns ++;
    }
  }
  /* validity check for references */
  for (i = 0; i < Np; i++) {
    if (photcode[i].type == PHOT_DEP) {
      Npri = photcodes[0].hashcode[photcode[i].equiv];
      if ((Npri >= Np) || (Npri < 0)) {
	fprintf (stderr, "reference for dependent photcode is not in photcodes\n");
	free (photcode);
	return (FALSE);
      }
      if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) {
	fprintf (stderr, "reference for dependent photcode is not a primary or secondary code\n");
	free (photcode);
	return (FALSE);
      }
    }
    if (photcode[i].type == PHOT_ALT) {
      Npri = photcodes[0].hashcode[photcode[i].code];
      if ((Npri >= Np) || (Npri < 0)) {
	fprintf (stderr, "reference for alternate photcode is not in photcodes\n");
	free (photcode);
	return (FALSE);
      }
      if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) {
	fprintf (stderr, "reference for alternate photcode is not a primary or secondary code\n");
	free (photcode);
	return (FALSE);
      }
    }
  }
  photcodes[0].code = photcode;
  photcodes[0].Ncode = Np;

  return (TRUE);

}

void ParseColorTerms (char *terms, float *X, int *N) {

  int i;
  char *p;

  p = terms;

  for (i = 0; (p != NULL) && (i < NCTERMS); i++) {
    X[i] = atof (p);
    p = strchr (p, ',');
    if (p == (char *) NULL) continue;
    p ++;
  }
  *N = i;
}

/********** photcode lookups **********/

/* return photcode for given name */
PhotCode *GetPhotcodebyName (char *name) {
  
  int i;

  if (name == (char *) NULL ) return (NULL);
  
  for (i = 0; i < photcodes[0].Ncode; i++) {
    if (!strcmp (photcodes[0].code[i].name, name)) {
      return (&photcodes[0].code[i]);
    }
  }
  return (NULL);
}
/* return photcode.code for given name */
int GetPhotcodeCodebyName (char *name) {
  
  int i;
  
  if (name == (char *) NULL ) return (0);

  for (i = 0; i < photcodes[0].Ncode; i++) {
    if (!strcmp (photcodes[0].code[i].name, name)) {
      return (photcodes[0].code[i].code);
    }
  }
  return (0);
}
/* return equivalent photcode for given name */
PhotCode *GetPhotcodeEquivbyName (char *name) {
  
  int i, equiv;
  
  if (name == (char *) NULL ) return (NULL);

  for (i = 0; i < photcodes[0].Ncode; i++) {
    if (!strcmp (photcodes[0].code[i].name, name)) {
      if (photcodes[0].code[i].equiv == 0) return (NULL);
      equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv];
      if (equiv == -1) return (NULL);
      return (&photcodes[0].code[equiv]);
    }
  }
  return (NULL);
}
/* return equivalent photcode.code for given name */
int GetPhotcodeEquivCodebyName (char *name) {
  
  int i, equiv;
  
  if (name == (char *) NULL ) return (0);

  for (i = 0; i < photcodes[0].Ncode; i++) {
    if (!strcmp (photcodes[0].code[i].name, name)) {
      if (photcodes[0].code[i].equiv == 0) return (0);
      equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv];
      if (equiv == -1) return (0);
      return (photcodes[0].code[equiv].code);
    }
  }
  return (0);
}

/* return photcode for given code */
PhotCode *GetPhotcodebyCode (int code) {
  
  int entry;
  
  if (code < 0) return (NULL);
  if (code > 0x10000) return (NULL);

  entry = photcodes[0].hashcode[code];
  if (entry == -1) return (NULL);

  return (&photcodes[0].code[entry]);
}
/* return photcode.code for given code */
char *GetPhotcodeNamebyCode (int code) {
  
  int entry;
  
  if (code < 0) return (NULL);
  if (code > 0x10000) return (NULL);

  entry = photcodes[0].hashcode[code];
  if (entry == -1) return (NULL);

  return (photcodes[0].code[entry].name);
}
/* return equivalent photcode for given code */
PhotCode *GetPhotcodeEquivbyCode (int code) {
  
  int entry, equiv;
  
  if (code < 0) return (NULL);
  if (code > 0x10000) return (NULL);

  entry = photcodes[0].hashcode[code];
  if (entry == -1) return (NULL);

  if (photcodes[0].code[entry].equiv == 0) return (NULL);
  equiv = photcodes[0].hashcode[photcodes[0].code[entry].equiv];

  if (equiv == -1) return (NULL);
  return (&photcodes[0].code[equiv]);
}
/* return equivalent photcode.code for given code */
int GetPhotcodeEquivCodebyCode (int code) {
  
  int entry;
  
  if (code < 0) return (0);
  if (code > 0x10000) return (0);

  entry = photcodes[0].hashcode[code];
  if (entry == -1) return (0);
  return (photcodes[0].code[entry].equiv);
}

int GetPhotcodeNsec (int code) {
  
  int Nsec;
  
  if (code < 0) return (-1);
  if (code > 0x10000) return (-1);

  Nsec = photcodes[0].hashNsec[code];
  return (Nsec);
}

/* Nsec of 0 is PRI */
PhotCode *GetPhotcodebyNsec (int Nsec) {
  
  if (Nsec > photcodes[0].Nsecfilt) return (NULL);
  if (Nsec < 0) return (NULL);
  
  return (&photcodes[0].code[Nseclist[Nsec]]);
}

int GetPhotcodeNsecfilt () {
  return (photcodes[0].Nsecfilt);
}

/* ALLOCATE and return list of all photcodes
   with photcode.equiv == code */
int *GetPhotcodeEquivList (int code, int *nlist) {

  int i, Nlist;
  int *list;

  ALLOCATE (list, int, MAX (1, photcodes[0].Ncode));
  Nlist = 0;
  for (i = 0; i < photcodes[0].Ncode; i++) {
    if (photcodes[0].code[i].equiv != code) continue;
    list[Nlist] = photcodes[0].code[i].code;
    Nlist ++;
  }
  REALLOCATE (list, int, MAX (1, Nlist));

  *nlist = Nlist;
  return (list);
}

/**** conversion to INTERNAL vs TABLE types makes the iPhot versions irrelevant ****/
# if (0) 
/******** photometry conversions *********/
double PhotInst (Measure *measure) {

  short Mi;
  double M;

  Mi = iPhotInst (measure);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotCat (Measure *measure) {

  short Mi;
  double M;

  Mi = iPhotCat (measure);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotSys (Measure *measure, Average *average, SecFilt *secfilt) {

  short Mi;
  double M;

  Mi = iPhotSys (measure, average, secfilt);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotRel (Measure *measure, Average *average, SecFilt *secfilt) {

  short Mi;
  double M;

  Mi = iPhotRel (measure, average, secfilt);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {

  short Mi;
  double M;

  Mi = iPhotCal (thisone, average, secfilt, measure, code);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {

  short Mi;
  double M;

  Mi = iPhotRef (code, average, secfilt, measure);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) {

  short Mi;
  double M;

  Mi = iPhotAve (code, average, secfilt);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) {

  short Mi;
  double M;

  Mi = iPhotdM (code, average, secfilt);
  M = 0.001*Mi;
  return (M);
}

/****/
double PhotXm (PhotCode *code, Average *average, SecFilt *secfilt) {

  int Mi;
  double Xm;

  Mi = iPhotXm (code, average, secfilt);
  Xm = (Mi == NO_MAG) ? -1.0 : pow (10.0, 0.01*Mi);
  return (Xm);
}

# endif 

/******** internal photometry conversions (keeps values in short int millimags) *********/
float PhotInst (Measure *measure) {

  int Np;
  float M;

  Np = photcodes[0].hashcode[measure[0].source];
  if (Np == -1) return (NO_MAG_PS);

  if (photcodes[0].code[Np].type == PHOT_REF) {
    M = measure[0].M_PS;
    return (M);
  }

  M = measure[0].M_PS - measure[0].dt_PS - ZERO_POINT;
	  
  return (M);

}

float PhotCat (Measure *measure) {

  int Np;
  float Mcat;
  PhotCode *code;

  Np = photcodes[0].hashcode[measure[0].source];
  if (Np == -1) return (NO_MAG_PS);

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Mcat = measure[0].M_PS;
    return (Mcat);
  }
  code = &photcodes[0].code[Np];
  Mcat = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C;
  
  return (Mcat);
}

float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) {

  int i, Np;
  float Mcat, Mcol, Msys, mc, Mc;
  PhotCode *code;

  Np = photcodes[0].hashcode[measure[0].source];
  if (Np == -1) return (NO_MAG_PS);

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Msys = measure[0].M_PS;
    return (Msys);
  }
  code = &photcodes[0].code[Np];
  Mcat = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C;

  /* for DEP, color must be made of PRI/SEC */
  mc = iPhotColor (average, secfilt, NULL, code);
  if (mc == NO_MAG_PS) return (Mcat);
  mc = mc - F_PS*code[0].dX;

  Mc = mc;
  Mcol = 0;
  for (i = 0; i < code[0].Nc; i++) {
    Mcol += code[0].X[i]*Mc;
    Mc *= mc;
  }
  Msys = Mcat + Mcol;
  return (Msys);
}

float PhotRel (Measure *measure, Average *average, SecFilt *secfilt) {

  int i, Np;
  float Mcat, Mcol, Mrel, mc, Mc;
  PhotCode *code;

  Np = photcodes[0].hashcode[measure[0].source];
  if (Np == -1) return (NO_MAG_PS);

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Mcat = measure[0].M_PS;
    return (Mcat);
  }
  code = &photcodes[0].code[Np];
  Mrel = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C - measure[0].Mcal_PS;

  /* for DEP, color must be made of PRI/SEC */
  mc = iPhotColor (average, secfilt, NULL, code);
  if (mc == NO_MAG_PS) return (Mrel);
  mc = mc - F_PS*code[0].dX;

  Mc = mc;
  Mcol = 0;
  for (i = 0; i < code[0].Nc; i++) {
    Mcol += code[0].X[i]*Mc;
    Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
  }
  Mrel += Mcol;
  return (Mrel);
}

/* return calibrated magnitude from measure for given photcode */
float PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {

  int i, Np; 
  float Mcal, Mrel, Mcol, mc, Mc;

  /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */
  Np = photcodes[0].hashcode[thisone[0].source];
  if (Np == -1) {
    return (NO_MAG_PS);
  }

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Mrel = thisone[0].M_PS;
    return (Mrel);
  }
  if (code[0].code != photcodes[0].code[Np].equiv) {
    return (NO_MAG_PS);
  }

  Mcal = PhotRel (thisone, average, secfilt) + F_PS*code[0].C;

  mc = iPhotColor (average, secfilt, measure, code);
  if (mc == NO_MAG_PS) return (Mcal);
  mc = mc - F_PS*code[0].dX;

  Mc = mc;
  Mcol = 0;
  for (i = 0; i < code[0].Nc; i++) {
    Mcol += code[0].X[i]*Mc;
    Mc *= mc;
  }
  Mcal += Mcol;
  return (Mcal);
}

/* color term may not use DEP magnitude */
float iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {

  int i, Ns1, Ns2, Ns;
  float m1, m2, mc;
  PhotCode *color;

  m1 = m2 = NO_MAG_PS;

  if (measure == NULL) {
    Ns1 = photcodes[0].hashNsec[code[0].c1];
    Ns2 = photcodes[0].hashNsec[code[0].c2];
  
    m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M_PS;
    m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M_PS;
    mc = ((m1 == NO_MAG_PS) || (m2 == NO_MAG_PS)) ? NO_MAG_PS : (m1 - m2);
    return (mc);
  }

  /* find magnitude matching first color term */
  color = GetPhotcodebyCode (code[0].c1);
  if (color == NULL) return (NO_MAG_PS);
  if (color[0].type == PHOT_REF) {
    for (i = 0; (i < average[0].Nm) && (m1 == NO_MAG_PS); i++) {
      if (measure[i].source == color[0].code) {
	m1 = measure[i].M_PS;
      }
    }
  } else {
    Ns = photcodes[0].hashNsec[color[0].code];
    m1 = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  }	

  /* find magnitude matching second color term */
  color = GetPhotcodebyCode (code[0].c2);
  if (color == NULL) return (NO_MAG_PS);
  if (color[0].type == PHOT_REF) {
    for (i = 0; (i < average[0].Nm) && (m2 == NO_MAG_PS); i++) {
      if (measure[i].source == color[0].code) {
	m2 = measure[i].M_PS;
      }
    }
  } else {
    Ns = photcodes[0].hashNsec[color[0].code];
    m2 = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  }	
  mc = ((m1 == NO_MAG_PS) || (m2 == NO_MAG_PS)) ? NO_MAG_PS : (m1 - m2);
  return (mc);
}

/* return calibrated magnitude from average/secfilt for given photcode */
float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {

  int i, Ns;
  float Mave, Mref, Mcol, mc;
  double Mc;

  Ns = photcodes[0].hashNsec[code[0].code];
  Mave = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  Mref = Mave + F_PS*code[0].C;

  mc = iPhotColor (average, secfilt, measure, code);
  if (mc == NO_MAG_PS) return (Mref);
  mc = mc - F_PS*code[0].dX;

  Mc = mc;
  Mcol = 0;
  for (i = 0; i < code[0].Nc; i++) {
    Mcol += code[0].X[i]*Mc;
    Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
  }
  Mref += Mcol;
  return (Mref);
}

/***/
float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) {

  int Ns;
  float Mave;

  Ns = photcodes[0].hashNsec[code[0].code];
  Mave = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  return (Mave);
}

/*** note that this is NOT a wrapper around iPhotdM ***/
float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) {

  int Ns;
  float dM;

  Ns = photcodes[0].hashNsec[code[0].code];
  dM  = (Ns == -1) ? average[0].dM : secfilt[Ns].dM_PS;
  return (dM);
}

/*** note that this is NOT a wrapper around iPhotXm ***/
float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt) {

  int Ns;
  short Mi;
  float Xm;

  Ns = photcodes[0].hashNsec[code[0].code];
  Mi = (Ns == -1) ? average[0].Xm : secfilt[Ns].Xm;
  Xm = (Mi == NO_MAG_PS) ? -1.0 : pow (10.0, 0.01*Mi);
  return (Xm);
}

/* given a photcode pair c1 & c2, return the color of this star (NaN if not found) */
int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color) {

  int i, Ns;
  double M1, M2, dM;
  PhotCode *code;
  
  code = GetPhotcodebyCode (c1);
  if (code == NULL) return (FALSE);
  if (code[0].type == PHOT_REF) {
    for (i = 0; i < average[0].Nm; i++) {
      if (measure[i].source == c1) {
	M1 = measure[i].M_PS;
	goto filter1;
      }
    }	
    return (FALSE);
  } else {
    Ns = photcodes[0].hashNsec[code[0].code];
    M1 = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  }	

filter1:
  code = GetPhotcodebyCode (c2);
  if (code == NULL) return (FALSE);
  if (code[0].type == PHOT_REF) {
    for (i = 0; i < average[0].Nm; i++) {
      if (measure[i].source == c2) {
	M2 = measure[i].M_PS;
	goto filter2;
      }
    }	
    return (FALSE);
  } else {
    Ns = photcodes[0].hashNsec[code[0].code];
    M2 = (Ns == -1) ? average[0].M : secfilt[Ns].M_PS;
  }	
  
filter2:

  dM = M1 - M2;
  *color = dM;
  
  return (TRUE);
}

/* photcode table should have the following format: 

# code name     type  zero  airmass  offset  c1 c2  slope  <color>  primary
1    B        pri   24.0  0.15     -       -  -   -      -        -
2    B        pri   24.0  0.15     -       -  -   -      -        -
3    B1       sec   22.5  0.18     0.15    1  2   0.10   0.50     1
1000 USNO_B   ref   -     -        -       -  -   -      -        -

*/


/*
  Nc1 = photcodes[0].code[Np].c1;
  Ns1 = photcodes[0].hashNsec[Nc1];

  Nc2 = photcodes[0].code[Np].c2;
  Ns2 = photcodes[0].hashNsec[Nc2];

  Xlam = photcodes[0].code[Np].X[0];
  Klam = photcodes[0].code[Np].K;
  
  m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M;
  m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M;
*/
