# include <dvo.h>

/* we use a static variable to save the pre-loaded the photcodes.
   photcode conversion functions refer to this table for ref values */

/* We have three indexes:

   table[0].hashcode provides the table photcode sequence for a photcode value
   table[0].hashNsec provides the Nsec sequence for a photcode value
   table[0].codeNsec provides the photcode value for an Nsec sequence

   given photcode = table[0].code[i] and code = photcode[0].code

   hashcode[code] = i
   hashNsec[code] = Nsec or -1 if not a PRI/SEC code
   codeNsec[Nsec] = code
*/

static PhotCodeData *photcodes = NULL;

PhotCodeData *GetPhotcodeTable () {

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

/* the static ZERO_POINT is used by programs as an approximate nominal */
static double ZERO_POINT;

void SetZeroPoint (double ZP) {
  ZERO_POINT = ZP;
}

# define SCALE 0.001

/********** photcode lookup functions **********/

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

  if (name == 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 == 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 == 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 == 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) {
  
  int Ncode;
  int Nseq;

  if (Nsec > photcodes[0].Nsecfilt) return (NULL);
  if (Nsec < 0) return (NULL);
  
  Ncode = photcodes[0].codeNsec[Nsec];
  if (Ncode < 0) return (NULL);

  Nseq = photcodes[0].hashcode[Ncode];
  return (&photcodes[0].code[Nseq]);
}

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);
}

/******** photometry conversion functions *********/
float PhotInst (Measure *measure) {

  int Np;
  float M;

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

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

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

}

float PhotCat (Measure *measure) {

  int Np;
  float Mcat;
  PhotCode *code;

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

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Mcat = measure[0].M;
    return (Mcat);
  }
  code = &photcodes[0].code[Np];
  Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*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].photcode];
  if (Np == -1) return (NO_MAG);

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

  /* for DEP, color must be made of PRI/SEC */
  mc = PhotColorForCode (average, secfilt, NULL, code);
  if (mc == NO_MAG) return (Mcat);
  mc = mc - SCALE*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].photcode];
  if (Np == -1) return (NO_MAG);

  if (photcodes[0].code[Np].type == PHOT_REF) {
    Mcat = measure[0].M;
    return (Mcat);
  }
  code = &photcodes[0].code[Np];
  Mrel = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal;

  /* for DEP, color must be made of PRI/SEC */
  mc = PhotColorForCode (average, secfilt, NULL, code);
  if (mc == NO_MAG) return (Mrel);
  mc = mc - SCALE*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].photcode];
  if (Np == -1) {
    return (NO_MAG);
  }

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

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

  mc = PhotColorForCode (average, secfilt, measure, code);
  if (mc == NO_MAG) return (Mcal);
  mc = mc - SCALE*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 PhotColorForCode (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {

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

  m1 = m2 = NO_MAG;

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

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

  /* find magnitude matching second color term */
  color = GetPhotcodebyCode (code[0].c2);
  if (color == NULL) return (NO_MAG);
  if (color[0].type == PHOT_REF) {
    for (i = 0; (i < average[0].Nm) && (m2 == NO_MAG); i++) {
      if (measure[i].photcode == color[0].code) {
	m2 = measure[i].M;
      }
    }
  } else {
    Ns = photcodes[0].hashNsec[color[0].code];
    m2 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
  }	
  mc = ((m1 == NO_MAG) || (m2 == NO_MAG)) ? NO_MAG : (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) ? NO_MAG : secfilt[Ns].M;
  Mref = Mave + SCALE*code[0].C;

  mc = PhotColorForCode (average, secfilt, measure, code);
  if (mc == NO_MAG) return (Mref);
  mc = mc - SCALE*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) ? NO_MAG : secfilt[Ns].M;
  return (Mave);
}

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

  int Ns;
  float dM;

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

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

  int Ns;
  short Mi;
  float Xm;

  Ns = photcodes[0].hashNsec[code[0].code];
  Mi = (Ns == -1) ? NO_MAG : secfilt[Ns].Xm;
  Xm = (Mi == NO_MAG) ? -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].photcode == c1) {
	M1 = measure[i].M;
	goto filter1;
      }
    }	
    return (FALSE);
  } else {
    Ns = photcodes[0].hashNsec[code[0].code];
    M1 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
  }	

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].photcode == c2) {
	M2 = measure[i].M;
	goto filter2;
      }
    }	
    return (FALSE);
  } else {
    Ns = photcodes[0].hashNsec[code[0].code];
    M2 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
  }	
  
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;
*/
