/* mk4time.c
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Computes the time of a feature at a specified position angle,
 * based on the direction of the scan, and the scan duration.
 *
 * The direction and duration values are obtained from the FITS header(s).
 * The user need only supply the angle(s) of interest, 
 * and the image file name (FITS).
 *
 * Uses the "cfitsio" library to extract information from the FITS headers.
 * NOTE: The direction, scan time are obtained from the extension HDU.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Author: Alice Lecinski   HAO/NCAR   findtime.c [Mk3]
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Update: Andrew L. Stanger   HAO/NCAR   mk4time.c   21 December 2000
 * 26 Aug 2004: Added more comments. 
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 */

#include <stdio.h>
#include <ctype.h>
#include "math.h"
#include "fitsio.h"

#define FOREVER 1
#define DONEXT  -999
#define QUIT    0

char *file_name ;
int  Anglesflag, angleflag ;

int pmode = 1 ;

/* -------------------------------------------------------------------------- */

main (an, av)

int an ;
char *av [] ;
{

    float get_frac_offset () ;
    void prsum () ;
    void printerror () ;

    char *w ;		/* pointer to command line arguments */
    static char comment [80] ;
    static char date_obs [24] ;		/* Observation date. */
    static char time_obs [24] ;		/* Observation time. */

    int bar_dir ;			/* Barrel scan direction. */
					/* +1 forward [CCW], -1 reverse [CW]. */
    int doit ;
    int phour, pminute, psecond ;	/* Position time. */
    int bhour, bminute, bsecond ;	/* Begin scan time. */
    int dhour, dminute, dsecond ;	/* Delta (offset) time. */
    int byear, bmonth,  bday    ;	/* Begin date. */
    int pyear, pmonth,  pday    ;	/* Position date. */
    int status ;

    float ang_range = 360.0 ;		/* Range of angles in scan. */
    float delta_time ;			/* Delta time (offset from start). */
    float duration ;			/* Scan duration (seconds). */
    float frac_time ;			/* Fraction of time for position. */
    float Phase ;			/* Phase angle (offset from Geo north.*/
    float raw_tmin, raw_tmax ;		/* Raw min/max barrel angles. */
    float solar_p0 ;			/* P-angle. */

    double PA ;
    double thetamin , thetamax , ThetaRng , fabs() ;

    FILE *infid ;
    fitsfile *fptr ;

    /*** Process command line arguments. ***/

    if ( an <= 2 ) { prsum () ; exit (-4) ; }
    while (--an)
	{
	if (*(w = *++av) == '-')	/* Check for flag arguments */
	    switch (w [1])
	    {
	    case 'A':  Anglesflag++ ; break ;
	    case 'a':  angleflag++ ;
	               --an ; w = *++av ;
		       sscanf (w, "%lf", &PA) ; break ;
	    default:   fprintf (stderr," unknown flag %s\n",*av) ;
		       prsum () ; exit (-1) ;
	    }

	else		/* Get the image filename. */
	    {

	    if (*w == '?') { prsum () ; exit (-2) ; }
	    if ( angleflag + Anglesflag == 0 )
		{
		fprintf (stderr,"One of -a or -A  MUST be given\n") ;
	        prsum () ; exit (-3) ;
		}

	    file_name = w ;
	    status = 0;

	    /*** Open FITS file. ***/

	    if ( fits_open_file (&fptr, file_name, READONLY, &status) )
		printerror ( status );

	    /*** Primary FITS header. ***/

	    if (fits_read_key (fptr, TSTRING, "DATE-OBS", date_obs,
		comment, &status))  printerror (status) ;
	    sscanf (date_obs, "%04d-%02d-%02d", &byear, &bmonth, &bday) ;
	    pyear  = byear ;
	    pmonth = bmonth ;
	    pday   = bday ;

	    if (fits_read_key (fptr, TSTRING, "TIME-OBS", time_obs,
		comment, &status))  printerror (status) ;
	    sscanf (time_obs, "%2d:%2d:%2d", &bhour, &bminute, &bsecond) ;

	    if (fits_read_key (fptr, TFLOAT, "SOLAR_P0", &solar_p0,
		comment, &status))  printerror (status) ;

	    /*** Move to extension header. ***/

	    if (fits_movabs_hdu (fptr, 2, NULL, &status))  printerror (status) ;

	    if (fits_read_key (fptr, TINT, "BAR_DIR", &bar_dir,
		comment, &status))  printerror (status) ;

	    if (fits_read_key (fptr, TFLOAT, "DURATION", &duration,
		comment, &status))  printerror (status) ;

	    if (fits_read_key (fptr, TFLOAT, "RAW_TMIN", &raw_tmin,
		comment, &status))  printerror (status) ;

	    if (fits_read_key (fptr, TFLOAT, "RAW_TMAX", &raw_tmax,
		comment, &status))  printerror (status) ;

	    if (fits_read_key (fptr, TFLOAT, "PHASE", &Phase,
		comment, &status))  printerror (status) ;

	    duration = duration / 60.0 ;	/* Convert from Hertz to sec. */

	    if (bar_dir > 0)
	        ang_range = raw_tmax - raw_tmin ;
	    else
		ang_range = raw_tmin - raw_tmax ;

	    if (pmode)
		{
	        printf ("raw_tmin/raw_tmax: %7.2f  %7.2f\n",
		         raw_tmin, raw_tmax) ;
		
		printf ("ang_range: %7.2f  duration: %7.2f\n",
			 ang_range, duration) ;
		}

	    doit = FOREVER ;

	    if ( Anglesflag )
		fprintf (stdout, 
		"Enter %d as a position angle to get to the next image.\n", 
		DONEXT ) ;

	    /*** Multiple position angle loop. ***/

	    while ( doit )
		{
		if ( Anglesflag )
		    {
		    fprintf (stdout, 
		       "\nEnter in the position angle of your feature --> ") ;
		    fscanf(stdin, "%lf", &PA) ;
		    if ( PA == DONEXT ) break ;
		    }
		else doit = QUIT ;

		/*** Compute fractional time offset from beginning of scan. ***/

		frac_time = get_frac_offset ( bar_dir,  ang_range,
					      solar_p0, Phase,     PA ) ;

		delta_time = duration * frac_time ;
		dhour   = (int)  (delta_time / 3600.0 );
		dminute = (int) ((delta_time - dhour   * 3600) / 60.0)  ;
		dsecond = (int)  (delta_time - dminute *   60) ;

		psecond = bsecond + dsecond ;
		pminute = bminute + dminute ;
		phour   = bhour   + dhour ;

		/*** Adjust time fields for overflow. ***/

		if (psecond > 59)
		    {
		    pminute += 1 ;
		    psecond -= 60 ;
		    }
		if (pminute > 59)
		    {
		    phour += 1 ;
		    pminute -= 60 ;
		    }
		if (phour > 23)
		    {
		    pday += 1 ;
		    phour -= 24 ;
		    }

                fprintf (stdout,
			 "%s  %04d-%02d-%02d  %02d:%02d:%02d @ %.2lf degrees\n",
		         file_name,
			 pyear, pmonth, pday, phour, pminute, psecond, PA) ;

		if (pmode)
		    fprintf (stdout,
			     "dir:%2d  delta_time:%.2lf  frac_time:%.2lf\n",
		              bar_dir, delta_time, frac_time) ;

		if (pmode)
		    printf ("\n") ;
	        }

	    /*** Close FITS file. ***/

	    if ( fits_close_file (fptr, &status) )  printerror (status) ;
	    }
        }
   exit (0) ;
}

/* -------------------------------------------------------------------------- */

/* get_frac_offset ()
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Compute fractional time offset from beginning of scan.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 */

float
get_frac_offset (dir, range, pangle, phase, position)

int dir ;			/* Scan direction: +1 CCW, -1 CW. */
float range ; 			/* Barrel angle range [degrees]. */
float pangle ;			/* Solar P-angle. */
float phase ;			/* Mk4 detector position angle. */
float position ;		/* Position Angle [CCW from solar north]. */
{

    float frac_offset ;
    float ref ;

    static char scan_dir [4];

    if (dir > 0) sprintf (scan_dir, "CCW");
    if (dir < 0) sprintf (scan_dir, "CW");

    /*** Forward scan direction: CCW [counter-clockwise]. ***/

    if ( dir > 0 )
	{
	ref = position - 180 + pangle + phase ;
        }

    /*** Reverse scan direction:  CW [clockwise]. ***/

    if ( dir < 0 )
	{
	ref = 360 - ( position - 180 + pangle + phase  ) ;
        }

    /*** Adjust angle so that it is within the range: 0 to 360. ***/

    if ( ref <   0 ) ref += 360 ;
    if ( ref > 360 ) ref -= 360 ;

    /*** Compute fractional offset. ***/

    frac_offset = (ref / range) ;

    if (pmode)
	{
	printf ("dir: %2d  %s  range: %7.2f  pangle: %7.2f  phase: %7.2f\n",
		dir, scan_dir, range, pangle, phase);
        printf ("position: %7.2f  ref: %7.2f  frac_offset: %7.2f\n",
		position, ref, frac_offset);
        }

/*
 *    printf ("ref:%.2lf  position:%.2lf  pangle:%.2lf  phase:%.2lf\n",
 *	     ref, position, pangle, phase) ;
 */

/*
 *    printf ("ref: %7.2f  range: %7.2f   frac_offset: %8.3f\n",
 *	      ref, range, frac_offset) ;
 */

    return (frac_offset) ;
}

/* -------------------------------------------------------------------------- */

/* prsum ()
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Print usage summary.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 */

void
prsum ()
{
    static char help [] = "\
mk4time [-a position_angle] [-A]  imagename1 imagename2 imagename3 \n\
\n\
Example:              findtime -a 270 19_20.srq \n\
Example:              findtime -a 270 19_07.srq 19_20.srq \n\
Example:              findtime -A     19_07.srq 19_20.srq \n\
\n\
This program attempts to find the exact time for a feature \n\
at the specified position angle in a Mk4 image.  \n\
\n\
Flags	    Description\n\
-----	    -----------\n\
\n\
-[a angle]  The position angle of interest.\n\
-[A]        Program will prompt the user for multiple position angles.\n\
" ;

    fputs (help, stderr) ;
}

/* -------------------------------------------------------------------------- */

void 
printerror ( int status )
{
    /*****************************************************/
    /* Print out cfitsio error messages and exit program */
    /*****************************************************/


    if (status)
        {
        fits_report_error (stderr, status); /* print error report */

        exit ( status );    /* terminate the program, returning error status */
        }

    return;
}

/* -------------------------------------------------------------------------- */