# include "astro.h"

int specpairfit (int argc, char **argv) {
  
  int i;
  Vector *flux1, *flux2, *dflux1, *dflux2, *window;

  if (argc != 7) goto usage;

  if ((flux1  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
  if ((dflux1 = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
  if ((flux2  = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
  if ((dflux2 = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) goto escape;
  if ((window = SelectVector (argv[5], OLDVECTOR, TRUE)) == NULL) goto escape;
  
  // XXX enforce matching lengths on the 6 vectors

  int mask = atoi (argv[6]);
  CastVector (window, OPIHI_INT);

  // minimize (flux1 - flux2*alpha) in window defined by mask
  // note that the mask is a SELECTION mask not an EXCLUSION mask

  double F12 = 0.0; 
  double F22 = 0.0;
  for (i = 0; i < flux1->Nelements; i++) {
    if ((mask & window->elements.Int[i]) == 0) continue;
    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
    F12 += flux1->elements.Flt[i] * flux2->elements.Flt[i] * weight;
    F22 += flux2->elements.Flt[i] * flux2->elements.Flt[i] * weight;
  }

  double Ao = F12 / F22;
  double dA = sqrt(1.0 / F22);

  int Ndof = -1; // 1 parameter fit
  double chisq = 0.0;
  for (i = 0; i < flux1->Nelements; i++) {
    if ((mask & window->elements.Int[i]) == 0) continue;
    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
    chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i]) * weight;
    Ndof ++;
  }

  double chisqNu = chisq / Ndof;

  // fprintf (stderr, "Ao: %f +/- %f, chisq: %f, chisq_nu : %f for %d dof\n", Ao, dA, chisq, chisqNu, Ndof);
  set_variable ("Ao", Ao);
  set_variable ("dA", dA);
  set_variable ("Xv", chisqNu);
  set_variable ("Nd", Ndof);
  return (TRUE);
  
 escape: 
  gprint (GP_ERR, "invalid vector\n");
  return (FALSE);
  
usage:
  gprint (GP_ERR, "USAGE: specpairfit (flux1) (dflux1) (flux2) (dflux2) (window) (mask) [options]\n");
  return FALSE;
}
