/*      FILE: convert.c
 *   PURPOSE: Convert instrumental to standard magnitudes
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001JUL04
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include "mx.h"
#include "wfpc2color.h"

int wfpc2color2_Convert_i1 (
  wfpc2color2_ParS *ParS
){
  char mxfunc[] = "wfpc2color2_Convert_i1";
  int status = 0;
  
  int 
    j,
    k,
    xxi,
    yyi,
    ccdnumi,
    one,
    two,
    id1i,
    id2i,
    n1,
    n2,
    match_not_found,
    color_ok,
    index[3];
 
  double 
    electrons_per_dn[3],
    gain_ratio[3],
    electrons_halfarcsec[3],
    dn[3],
    dnpersec[3],
    wfpc2_mag[3],
    imagerr[3],
    imag[3],
    mag[3],
    magerr[3],
    xd[3],
    yd[3],
    color,
    color2,
    color_err,
    exponent,
    expsec,
    gain,
    cte,
    synmag[3],
    zsyn,
    aimag,
    a,
    b,
    c,
    q,
    root,
    root1,
    root2,
    signb,
    bb,
    ac4,
    t1[3],
    t2[3],
    z[3];

  char
    namec[3][14+1];

  double
    gain_ratios[5] = { 0.0, 1.987, 2.003, 2.006, 1.955};

  /* use the following constant instead of the approximation 2.512 */
  const double magnitude_constant = pow(10.0,0.4);

  wfpc2color2_ImagS 
    *S1,
    *S2;

  S1 = ParS->ImagS[1];
  S2 = ParS->ImagS[2];

  n1 = ParS->instars[1];
  n2 = ParS->instars[2];

for (j=1; j<3; j++) {
  t1[j] = ParS->t1[j];
  t2[j] = ParS->t2[j];
  z[j] = ParS->z[j];
}


for (k=0; k<n1; k++) {

  one = k;
  two = 0;
  id1i = S1[one].idi;
  match_not_found = MX_TRUE;
  while (match_not_found&&(two<n2)) {
    id2i = S2[two].idi;
    if (id2i==id1i) {
      match_not_found = MX_FALSE;
    } else {
      two++;
    }
  }
  if (match_not_found) goto next_star;
  index[1] = one;
  index[2] = two;
  
  for (j=1; j<3; j++){

    xd[j] = ParS->ImagS[j][index[j]].xd;
    yd[j] = ParS->ImagS[j][index[j]].yd;
    imag[j] = ParS->ImagS[j][index[j]].magd;
    imagerr[j] = ParS->ImagS[j][index[j]].errd;
    magerr[j] = sqrt( imagerr[j]*imagerr[j] + (0.01*0.01) );  
    if (magerr[j]>ParS->mxerr[j]) goto next_star;

    /* set proper gain ratios */
    gain_ratio[j] = 1.0;
    if (ParS->wfpc2_gain[j] == WFPC2_GAIN_LOW) 
      gain_ratio[j] = gain_ratios[ParS->ccdnum[j]];

    /* convert to ELECTRONS in the standard half-arcec diameter aperture */
    exponent = 
      imag[j]             /* instrumental magnitude */
      + (ParS->apcor[j])  /* aperture correction to halfarcsec diameter */
      + (ParS->delta[j]); /* breathing correction */
    electrons_halfarcsec[j] = pow(magnitude_constant,-exponent);

    /* convert to DN */
    gain = ParS->gain[j];
    if (gain<=0.0) { status = 2;
      sprintf (MX.tmpmsg, "#   gain values must be greater than zero!\n");
      goto mx_error2;
    }
    dn[j] = electrons_halfarcsec[j]/gain;

    /* convert to DN_PER_SEC */
    expsec = ParS->expsec[j];
    if (expsec<=0.0) { status = 3;
      sprintf (MX.tmpmsg, "#   exposure times must be greater than zero!\n");
      goto mx_error2;
    }
    dnpersec[j] = dn[j]/expsec;

    /* correct for CTE in the Y direction*/
    cte = ParS->ycte[j];
    if (cte<0.0) { status = 4;
      sprintf (MX.tmpmsg, "#   Y_CTE corrections can not be negative!\n");
      goto mx_error2;
    }
    dnpersec[j] *= 1.0 + (cte*(yd[j]/800.0));

    /* correct for CTE in the X direction*/
    cte = ParS->xcte[j];
    if (cte<0.0) { status = 4;
      sprintf (MX.tmpmsg, "#   X_CTE corrections can not be negative!\n");
      goto mx_error2;
    }
    dnpersec[j] *= 1.0 + (cte*(xd[j]/800.0));

    /* convert to WFPC2 magnitudes */
    wfpc2_mag[j] = -2.5*log10(dnpersec[j]) + 2.5*log10(gain_ratio[j]);

    /* compute synthetic WFPC2 magnitudes */
    zsyn = ParS->zsyn[j];
    synmag[j] = wfpc2_mag[j] + zsyn;

    /* remove the reddening due to the instrumental filters */
    aimag = ParS->aimag[j];
    wfpc2_mag[j] -= aimag;
      
  } /* for (j=1; j<3; j++) */

  if (ParS->do_transform) {
    if ( ParS->SIGN == -1) {
      a = t2[2] - t2[1];
      b = t1[2] - t1[1] - 1.0;
      c = (z[2] - z[1]) + (wfpc2_mag[2]-wfpc2_mag[1]);
    } else { 
      a = t2[1] - t2[2];
      b = t1[1] - t1[2] - 1.0;
      c = (z[1] - z[2]) + (wfpc2_mag[1]-wfpc2_mag[2]);
    }
    bb = b*b;
    ac4 = a*c*4.0;
    if ( a == 0.0 ) { /* no quadratic color term */
      if ( b == 0.0 ) goto next_star; /* color is undefined! */
      color = -c/b;
    } else {
      if (bb>ac4) {
	signb = 1.0;
	if (b<0.0) signb = -1.0;
	q = -(0.5)*(b + signb*sqrt(bb-ac4));
	root1 = q/a;
	root2 = c/q;
	root = root1;
	if (fabs(root1)>fabs(root2)) root = root2;
	color = root;
      } else {
	goto next_star;
      }
    }
    color2 = color*color;
    for (j=1; j<3; ++j) 
      mag[j] = wfpc2_mag[j] + (t1[j]*color) + (t2[j]*color2) + (z[j]);
    /* Add back reddening due to the standard filters */
    mag[1] += ParS->asmag[1];
    mag[2] += ParS->asmag[2]; 
  } else {
    mag[1] = synmag[1];
    mag[2] = synmag[2];
    ParS->SIGN = 1;
  }
  color = (mag[1] - mag[2])*ParS->SIGN;
  color_ok =
       (color>=ParS->cmin[1]) 
    && (color<=ParS->cmax[1]) 
    && (color>=ParS->cmin[2]) 
    && (color<=ParS->cmax[2])
  ;
#ifdef DEBUG
  fprintf(stderr,"%7.3f >= %7.3f <= %7.3f >= %7.3f <= %7.3f : %d\n",
    color,
    ParS->cmin[1],
    ParS->cmax[1],
    ParS->cmin[2],
    ParS->cmax[2],
    color_ok
  );
#endif
  if (!color_ok) goto next_star;
  color_err = sqrt( (magerr[1]*magerr[1]) + (magerr[2]*magerr[2]) );
  
  if (ParS->SIGN!=-1) {
    printf( "%8.3f%7.3f", mag[1], magerr[1]);
    printf( "%8.3f%7.3f", mag[2], magerr[2]);
  } else { /* switch columns */
    printf( "%8.3f%7.3f", mag[2], magerr[2]);
    printf( "%8.3f%7.3f", mag[1], magerr[1]);
  }
  printf( "%8.3f%7.3f ", color, color_err);
  printf( "%8.3f%8.3f", synmag[1], synmag[2]);

  printf( "%7d ", id1i);

  /*
   * =====================================================================
   */

  for (j=1; j<3; j++){
    ccdnumi = ParS->ccdnum[j];
    xxi = (ParS->ImagS[j][index[j]].xd + 0.05)*10.0;
    yyi = (ParS->ImagS[j][index[j]].yd + 0.05)*10.0;
    sprintf( namec[j], "%02d%06d%06d", ccdnumi, yyi, xxi);
  }
  printf( "%s %s", namec[1], namec[2]);

  for (j=1; j<3; j++) {
    printf( 
      "   %1d%8.2f%8.2f%8.3f%7.3f", 
      j,
      ParS->ImagS[j][index[j]].xd,
      ParS->ImagS[j][index[j]].yd,
      ParS->ImagS[j][index[j]].magd,
      ParS->ImagS[j][index[j]].errd
    );
  }

  printf( " \"%s\"", ParS->field);

  printf( "\n");

next_star:
;
}

  printf ("#    (1)");
  printf ("    (2)");
  printf ("     (3)");
  printf ("    (4)");
  printf ("     (5)");
  printf ("    (6)");
  printf ("      (7)");
  printf ("     (8)");
  printf ("    (9)");
  printf ("           (10)");
  printf ("           (11)");
  printf (" (12)");
  printf ("   (13)");
  printf ("    (14)");
  printf ("    (15)");
  printf ("   (16)");
  printf (" (17)");
  printf ("   (18)");
  printf ("    (19)");
  printf ("    (20)");
  printf ("   (21)");
  printf (" (22)");
  printf ("\n");
  printf ("#\n");
  printf ("# column description\n");
  if (ParS->SIGN!=-1) {
    printf ("#    (1) SMAG1 magnitude: \"%s\"\n",ParS->SMAG[1]); 
    printf ("#    (2) rms error of SMAG1 [mag]\n"); 
    printf ("#    (3) SMAG2 magnitude: \"%s\"\n",ParS->SMAG[2]); 
    printf ("#    (4) rms error of SMAG2 [mag]\n"); 
  } else { /* switch columns */
    printf ("#    (1) SMAG2 magnitude: \"%s\"\n",ParS->SMAG[2]); 
    printf ("#    (2) rms error of SMAG2 [mag]\n"); 
    printf ("#    (3) SMAG1 magnitude: \"%s\"\n",ParS->SMAG[1]); 
    printf ("#    (4) rms error of SMAG1 [mag]\n"); 
  }
  printf ("#    (5) SCOL: \"%s\" [mag]\n",ParS->SCOL); 
  printf ("#    (6) rms error of SCOL [mag]\n"); 
  printf ("#    (7) *synthetic* magnitude (IMAG1=\"%s\")\n",ParS->imag[1] );
  printf ("#    (8) *synthetic* magnitude (IMAG2=\"%s\")\n",ParS->imag[2] );
  printf ("#    (9) Integer identifier\n"); 
  printf ("#   (10) 1st image identifier [format: ChipnumYYY.YXXX.X]\n"); 
  printf ("#   (11) 2nd image identifier [format: ChipnumYYY.YXXX.X]\n"); 
  printf ("#   (12) FILE1 is \"%s\"\n",ParS->in1);
  printf ("#   (13) FILE1: X position\n"); 
  printf ("#   (14) FILE1: Y position\n"); 
  printf ("#   (15) FILE1: *instrumental* magnitude (IMAG1=\"%s\") \n",ParS->imag[1] ); 
  printf ("#   (16) FILE1: rms error of IMAG1\n"); 
  printf ("#   (17) FILE2 is \"%s\"\n",ParS->in2);
  printf ("#   (18) FILE2: X position\n"); 
  printf ("#   (19) FILE2: Y position\n"); 
  printf ("#   (20) FILE2: *instrumental* magnitude (IMAG2=\"%s\") \n",ParS->imag[2] ); 
  printf ("#   (21) FILE2: rms error of IMAG2\n"); 
  printf ("#   (22) field name\n"); 

ok:
  status = 0;
  goto bye;

mx_error1:
  mxp_tmpmsg_init_f0();
mx_error2:
  mxp_errmsg_append_f3 (mxfunc, status, MX.tmpmsg);
  mxp_tmpmsg_init_f0();
  status = 1;
  goto bye;

bye:
  return(status);

}
/* end-of-file */