/* pointsharc -- Suggest (FAZO, FZAO) for science target based on
                 pointing model and observations of pointing reference */

/* C. Darren Dowell */

/* version 1.0 -- 2006 Mar 3 */
/* version 1.1 -- 2007 Jan 9 -- quantify pointing drift */
/* version 1.2 -- 2007 Mar 23 -- fix bug in "before&after" result to stderr */

/* To compile under Linux/UNIX:  cc -o pointsharc pointsharc.c -lm */

#include <stdio.h>
#include <math.h>

#define  MAX_SCAN  100
#define  MAX_SCAN_NO  999999
#define  MAX_LEN  1000
#define  MINFWHM  5.0
#define  MAXFWHM  40.0
#define  MINVOLTS  0.03
#define  MAXVOLTS  4000.0

main(argc, argv)

int  argc;
char  *argv[];

{
  FILE  *modelfile, *calfile;
  int  tscan, nrscan, arg, rscan[MAX_SCAN], scan, i, j, n1, n2, n3;
  float  tmfazo, tmfzao, fazo, fzao;
  float  rmfazo[MAX_SCAN], rmfzao[MAX_SCAN], rcfazo[MAX_SCAN], rcfzao[MAX_SCAN];
  float  list1[MAX_SCAN][2], list2[MAX_SCAN][2], list3[MAX_SCAN][2];
  float  temp, dfazo1, dfzao1, dfazo2, dfzao2, dfazo3, dfzao3;
  float  dfazo, dfzao;
  char  line[MAX_LEN];

  /* READ COMMAND LINE ARGUMENTS */
  if (argc < 5)
    usage();
  if (sscanf(argv[3], "%d", &tscan) != 1)
    usage();
  if ((tscan < 0) || (tscan > MAX_SCAN_NO))
    printabort("Target scan number out of range");
  tmfazo = 1.0e30;
  tmfzao = 1.0e30;
  nrscan = 0;
  arg = 4;
  while (arg < argc) {
    if (sscanf(argv[arg], "%d", &(rscan[nrscan])) != 1)
      usage();
    if ((rscan[nrscan] < 0) || (rscan[nrscan] > MAX_SCAN_NO))
      printabort("Reference scan number out of range");
    rmfazo[nrscan] = 1.0e30;
    rmfzao[nrscan] = 1.0e30;
    rcfazo[nrscan] = 1.0e30;
    rcfzao[nrscan] = 1.0e30;
    nrscan++;
    arg++;
  }

  /* READ MODEL FILE */
  if ((modelfile = fopen(argv[1], "r")) == NULL)
    printabort("Couldn't open model file");
  while (fgets(line, MAX_LEN, modelfile) != NULL) {
    if (sscanf(line,
               "%d %*d %*d %*d %*f %*s %*f %*f %*f %*f %*f %*f %*s %f %f",
               &scan, &fazo, &fzao) == 3) {
      if (scan == tscan) {
        tmfazo = fazo;
        tmfzao = fzao;
      }
      for (i = 0; i < nrscan; i++) {
        if (scan == rscan[i]) {
          rmfazo[i] = fazo;
          rmfzao[i] = fzao;
        }
      }
    }
  }
  fclose(modelfile);

  /* READ CALIBRATOR RESULT FILE */
  if ((calfile = fopen(argv[2], "r")) == NULL)
    printabort("Couldn't open calibration reference file");
  while (fgets(line, MAX_LEN, calfile) != NULL) {
    if (sscanf(line,
               "%d %*d %*d %*d %*f %*f %*f %*f %f %f",
               &scan, &fazo, &fzao) == 3) {
      for (i = 0; i < nrscan; i++) {
        if (scan == rscan[i]) {
          rcfazo[i] = fazo;
          rcfzao[i] = fzao;
        }
      }
    }
  }
  fclose(calfile);
  
  /* DO THE CALCULATION */
  if (tmfazo > 1.0e29)
    printabort("No target model info.");
  n1 = n2 = n3 = 0;
  for (i = 0; i < nrscan; i++) {
    if (rmfazo[i] > 1.0e29) {
      fprintf(stderr, "WARNING:  No reference model info. for scan %d\n",
              rscan[i]);
    }
    else {
      if (rcfazo[i] > 1.0e29) {
        fprintf(stderr, "WARNING:  No reference cal. info. for scan %d\n",
                rscan[i]);
      }
      else {
        list1[n1][0] = rcfazo[i] - rmfazo[i];     /* "all" list */
        list1[n1][1] = rcfzao[i] - rmfzao[i];
        n1++;
        if (rscan[i] < tscan) {
          list2[n2][0] = rcfazo[i] - rmfazo[i];     /* "before" list */
          list2[n2][1] = rcfzao[i] - rmfzao[i];
          n2++;
        }
        if (rscan[i] > tscan) {
          list3[n3][0] = rcfazo[i] - rmfazo[i];     /* "after" list */
          list3[n3][1] = rcfzao[i] - rmfzao[i];
          n3++;
        }
      }
    }
  }
  for (i = 0; i < (n1-1); i++) {
    for (j = i+1; j < n1; j++) {
      temp = list1[i][0];
      if (list1[j][0] < temp) {
        list1[i][0] = list1[j][0];
        list1[j][0] = temp;
      }
    }
  }
  for (i = 0; i < (n1-1); i++) {
    for (j = i+1; j < n1; j++) {
      temp = list1[i][1];
      if (list1[j][1] < temp) {
        list1[i][1] = list1[j][1];
        list1[j][1] = temp;
      }
    }
  }
  for (i = 0; i < (n2-1); i++) {
    for (j = i+1; j < n2; j++) {
      temp = list2[i][0];
      if (list2[j][0] < temp) {
        list2[i][0] = list2[j][0];
        list2[j][0] = temp;
      }
    }
  }
  for (i = 0; i < (n2-1); i++) {
    for (j = i+1; j < n2; j++) {
      temp = list2[i][1];
      if (list2[j][1] < temp) {
        list2[i][1] = list2[j][1];
        list2[j][1] = temp;
      }
    }
  }
  for (i = 0; i < (n3-1); i++) {
    for (j = i+1; j < n3; j++) {
      temp = list3[i][0];
      if (list3[j][0] < temp) {
        list3[i][0] = list3[j][0];
        list3[j][0] = temp;
      }
    }
  }
  for (i = 0; i < (n3-1); i++) {
    for (j = i+1; j < n3; j++) {
      temp = list3[i][1];
      if (list3[j][1] < temp) {
        list3[i][1] = list3[j][1];
        list3[j][1] = temp;
      }
    }
  }
  dfazo1 = 0.0;
  dfzao1 = 0.0;
  if (n1 == 0) {
    fprintf(stderr, "WARNING:  No reference cal. information at all\n");
  }
  else {
    dfazo1 = (list1[n1/2][0] + list1[(n1-1)/2][0])/2.0;
    dfzao1 = (list1[n1/2][1] + list1[(n1-1)/2][1])/2.0;
  }
  dfazo2 = 0.0;
  dfzao2 = 0.0;
  if (n2 > 0) {
    dfazo2 = (list2[n2/2][0] + list2[(n2-1)/2][0])/2.0;
    dfzao2 = (list2[n2/2][1] + list2[(n2-1)/2][1])/2.0;
  }
  dfazo3 = 0.0;
  dfzao3 = 0.0;
  if (n3 > 0) {
    dfazo3 = (list3[n3/2][0] + list3[(n3-1)/2][0])/2.0;
    dfzao3 = (list3[n3/2][1] + list3[(n3-1)/2][1])/2.0;
  }
  dfazo = 0.0;
  dfzao = 0.0;
  if (n2 > 0) {
    if (n3 > 0) {
      dfazo = (dfazo2 + dfazo3)/2.0;
      dfzao = (dfzao2 + dfzao3)/2.0;
    }
    else {
      dfazo = dfazo2;
      dfzao = dfzao2;
    }
  }
  else {
    if (n3 > 0) {
      dfazo = dfazo3;
      dfzao = dfzao3;
    }
  }

  fprintf(stderr, "static pointing model:  (%.1f, %.1f)\n", tmfazo, tmfzao);
  fprintf(stderr, "'all' method:  dynamic residual = (%.1f, %.1f)\n",
          dfazo1, dfzao1);
  fprintf(stderr, "               suggested pointing = (%.1f, %.1f)\n",
          tmfazo + dfazo1, tmfzao + dfzao1);
  fprintf(stderr, "'before&after' method:  dynamic residual = (%.1f, %.1f)\n",
          dfazo, dfzao);
  fprintf(stderr, "                        suggested pointing = (%.1f, %.1f)\n",
          tmfazo + dfazo, tmfzao + dfzao);
  printf("%d  %.1f %.1f  %.1f %.1f", tscan, tmfazo +  dfazo1,
         tmfzao + dfzao1, tmfazo + dfazo, tmfzao + dfzao);
  if ((n2 > 0) && (n3 > 0)) {
    printf("  %4.1f %4.1f", dfazo3-dfazo2, dfzao3-dfzao2);
    fprintf(stderr, "pointing drift:  (%.1f, %.1f)\n", dfazo3-dfazo2, dfzao3-dfzao2);
  }
  printf("\n");
}

/****************************************************/

printabort(char  *message)

{
  fprintf(stderr, "ERROR:  %s\n", message);
  exit (1);
}

/****************************************************/

usage()

{
  fprintf(stderr, "Usage:  pointsharc modelfile pointreffile target_scan_# ref_scan_# {ref_scan_# ...}\n");
  exit (1);
}


