// sharctau.c
// Oct 23, 2003: Colin Borys
//
// This program reads the header information from a sharc2 fits file,
// as well as polynomial fits to the CSO tau made by Jon Bird, and
// computes the "fitted" tau at the start of the sharc2 observation.
//
//
//
// options are:
//    -v : verbose mode.  Reports basic header information.  Otherwise
//         program just spits out the fit tau.

//gcc -o sharctau sharctau.c -lcfitsio -lm


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/types.h>
#include <fcntl.h>

#include "fitsio.h"
#define MAX_STR_LEN 300

void usage();
void report_error( char *msg);
void report_fits_error( char *msg, int status);
void report_fits_header_error( char *fname, char *header,int status);
int getline(FILE* f,char *c);

int main(int argc, char *argv[]) {

  char sharcfn[MAX_STR_LEN]="sharc2-000000.fits";
  char fitsfn[MAX_STR_LEN]="fitfile.txt";
  char l[MAX_STR_LEN];   // a line of text
  FILE *fits;            // file pointer to fits file
  struct stat fbuffer;

  fitsfile *fp;                    // file pointer
  int      status;                 // status from fits command
  int i;                           // a counter
  int   isJSHARC;
  float TAU225GH;
  char  DATE[80];
  char  CREATOR[80];
 
  int year,month,day,hour,minute;

  float fut;                       // fractional UT
  int order;                       // order of polynomial to fit
  float coeff[20];                 // coefficients
  char utdate[MAX_STR_LEN];        // ut date in string format
  float utstart,utstop;            // range of valid UT times (fit range)
  int foundline;                   // found the line in the fits file?
  int done;                        // done reading the file?
  float fittau;                    // value of the fit tau;

  // set defaults to options
  int c;                           // for parsing options
  int verbose = 0;                 // turn verbose output off by default
  int nargrequired;                // number of arguments required.
  

  nargrequired=3;  // 3 arguments required (command + 2 filenames)
  while ((c=getopt(argc,argv,"v")) != EOF) {
    switch(c) {
    case 'v':
      verbose = 1;
      nargrequired++;
      break;
    default:
      fprintf(stderr,"Unrecognized option\n");
      usage();
    }
  }

  if (argc <nargrequired) usage();
  strncpy(sharcfn,argv[optind],MAX_STR_LEN);
  strncpy(fitsfn,argv[optind+1],MAX_STR_LEN);

  ////////////////////////////////////////////////////////
  // PART 1: Parse the data file

  // test to see if it is a "valid" file
  status=0; 
  status=stat(sharcfn,&fbuffer);
  if (status!=0) {  // some problem opening file!
    perror(argv[0]);
    exit(1);
  }
  if (fbuffer.st_size==0) {  // 0 length file!
    fprintf(stderr,"%s has a filelength of zero!\n",fitsfn);
    exit(1);
  }
 

  // open the file
  if (fits_open_file(&fp, sharcfn, READONLY, &status))
    report_fits_error(sharcfn,status);
     
  // IRC or JSHARC took the data?
  if (fits_read_key(fp,TSTRING,"CREATOR",   CREATOR, NULL, &status))
    report_fits_header_error(fitsfn,"CREATOR",status);
  if (strncmp(CREATOR, "JSharc", 6) == 0) isJSHARC=1; 
  else isJSHARC=0;


  if (fits_read_key(fp, TSTRING,"DATE",      DATE, NULL, &status))
    report_fits_header_error(fitsfn,"DATE",status);

  // need to do a date fix for JSharc created files.
  if (isJSHARC) {
    DATE[6]++;
    if (DATE[6] > '9') {
      DATE[6] = '0';
      DATE[5] ++;
    }
  }
  // get UT time and date
  if (fits_str2time(DATE,&year,&month,&day,&hour,&minute,NULL,&status))
    report_fits_header_error(fitsfn,"Can't convert Date",status);
  // get the reported 225GHz tau
  if (fits_read_key(fp, TFLOAT, "TAU225GH",&TAU225GH, NULL, &status))
    report_fits_header_error(fitsfn, "TAU225GH",status);

  // compute UT in units oİf fraction of a day.
  fut=(hour + minute/60.0)/24.0 ;

  //make the ut date string so we can search the fit file.
  snprintf(utdate,MAX_STR_LEN,"%4i%02i%02i",year,month,day);
  fits_close_file(fp,&status); 



  ////////////////////////////////////////////////////////
  // PART 2: Parse the fit file
  if((fits=fopen(fitsfn,"r"))==NULL) {
    perror(fitsfn);
    exit(1);
  }

  foundline=0;
  done=0;
  while (!done) {
    if (getline(fits,l)==0) done=1;
    else 
      if(strncmp(l,utdate,8)==0) {
	foundline=1;
	strtok(l," \t");
	sscanf(strtok(NULL," \t"),"%f",&utstart);
	sscanf(strtok(NULL," \t"),"%f",&utstop);
	sscanf(strtok(NULL," \t"),"%i",&order);
	for (i=0;i<=order;i++)
	  sscanf(strtok(NULL," \t"),"%f",&coeff[i]);
      }
  }
  fclose(fits);

  if (!foundline) 
    report_error("Error: Could not fine appropriate tau fit entry\n");


  //////////////////////////////////////
  // PART 3: compute fit value
  //
  if ((fut<utstart) || (fut>utstop)) 
    fprintf(stderr,"Warning: UT time outside range of fit function!\n");
  
  fittau=0;
  for (i=0;i<=order;i++) 
    fittau+=coeff[i]*pow(fut,i);

  if (verbose) {
    printf("# DD/MM/YYYY HH:MM UT/24 TAU225 FITTAU\n");
    printf("# %02i/%02i/%04i %02i:%02i %4.2f  %5.3f %5.3f\n",day,month,year, hour,minute,fut,TAU225GH, fittau);
  } else 
    printf("%5.3f\n",fittau);
  

  
  return(1);
}


void usage() {
  report_error("Usage:\n sharctau [-v] datafile taufile\n");
} 

void report_error( char *msg){
  fprintf(stderr,"%s ",msg);
  exit(1);
}

void report_fits_error( char *msg, int status){
  fprintf(stderr,"%s ",msg);
  fits_report_error( stderr, status);
  exit( status);
}

void report_fits_header_error( char *fname, char *header,int status){
  fprintf(stderr,"%s:%s",fname,header);
  fits_report_error( stderr, status);
  exit( status);
}

int getline(FILE* f,char *c) {
  // reads in a line from a file, returns the
  // number of characters, and null terminates 
  // the string 
  int i=0;          // number of characters read
  char a=0;
  while ((a!='\n') & (a!=EOF)) {
    a=fgetc(f);
    if (a!=EOF) {
      c[i]=a;
      i++;
    }
  }
  c[i]='\0';   // null terminate the string 
  return i;
}
