/*******************************************************************************
 * This program is a total hack of several programs.                           *
 * It computes the relative filling factor for different tracks.  It assumes   *
 *   that the track will be at roughly the ideal "best" az for that track.     *
 * This code makes use of the WCS package of fits readers and it must be       *
 *   installed to work.                                                        *
 *  Usage: checkfill fitsfilename                                              *
 *                                                                             *
 ******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <time.h>

#define R2D 57.29577958      /* degrees per radian */
#define P 0.493336           /* cos alt * cos latitude */
#define Q 0.417863           /* sin alt * sin latitude */
#define DXMAX 1800.          /* tracker limits in mm */
#define DXMIN -1800.
#define DYMAX 1900.
#define DYMIN -1900.
#define FBMAX 0.1484
#define FBMIN -0.1484
#define FTMAX 0.1484
#define FTMIN -0.1484
#define FPMAX 0.1484
#define FPMIN -0.1484
#define FS 12528         /* Fs in mm */
#define PI                 3.14159
#define rad                0.017453293
#define telangle          55.0
#define deltatelangle      8.485
#define latitude          30.6716667
#define longitude        104.0216667
#define beg(a,b)          (sin((telangle+b)*rad)-sin(latitude*rad)*sin(a*rad))/(cos(latitude*rad)*cos(a*rad))
#define end(a,b)          (sin((telangle-b)*rad)-sin(latitude*rad)*sin(a*rad))/(cos(latitude*rad)*cos(a*rad))
#define UT(a,b)           (a-100.46-0.9856474*b)/15.0
#define az(a)           acos((sin(a*rad)-sin(latitude*rad)*sin(telangle*rad))/(cos(latitude*rad)*cos(telangle*rad)))/rad

main(argc, argv)
     int argc;
     char *argv[];
{
	double r, d, x, w, z, beta, theta, phi, fs, c10, c20, c30, c40, c50, h;
	double ha, pa, hc, ha1, ha2, tde;
        double newstart, newend, bmin, xstrt, ystrt, obsaz;
        double betastart;
        double secperdeg, fill, fillsum;
	int inwindow, tstart, tend;
 int i, j, k, y, hr, min, deg;
 int tracks;
 double DECdeg, DECmin, DECsec;
 double DEC2, sec, asec, az;
 double tens, units, length, tmidpt, exp, ttime;
 double ha_start, ha_end;
 char in[75], temp[50], file[20];
 char DEC[25];
 char expo[10], xstart[20], ystart[20], actualaz[10];
 FILE *f1;

 if (argc != 2) {
    printf ("Usage: checkfill fitsfilename\n");
    exit ( 1 );
 }

 sscanf(argv[1],"%s%*c",&file);

 sprintf(temp,"gethead %s DEC STRUCTAZ X_STRT Y_STRT EXPTIME > checkfilltempfile",file );
 system(temp);
 f1=fopen("checkfilltempfile","r");
 fgets(in, 75, f1);
 fclose(f1);
 sscanf(in,"%s %s %s %s %s*c",&DEC,&actualaz,&xstart,&ystart,&expo);
 system("rm -f checkfilltempfile");

 for (i=0;i<strlen(actualaz);i++){
   if (actualaz[i]!='0' && actualaz[i]!='1' && actualaz[i]!='2' && actualaz[i]!='3' && actualaz[i]!='4' && actualaz[i]!='5' && actualaz[i]!='6' && actualaz[i]!='7' && actualaz[i]!='8' && actualaz[i]!='9' && actualaz[i]!='.') {
      printf("A data input ERROR has occured in the Azimuth; please go back.\n");
      exit(2);
   }
 }
 for (i=0;i<strlen(xstart);i++){
   if (xstart[i]!='0' && xstart[i]!='1' && xstart[i]!='2' && xstart[i]!='3' && xstart[i]!='4' && xstart[i]!='5' && xstart[i]!='6' && xstart[i]!='7' && xstart[i]!='8' && xstart[i]!='9' && xstart[i]!='.' && xstart[i]!='E' && xstart[i]!='-' && xstart[i]!='+') {
      printf(" A data input ERROR has occured in the X Start tracker position; please go back.\n");
      exit(2);
   }
 }
 for (i=0;i<strlen(ystart);i++){
   if (ystart[i]!='0' && ystart[i]!='1' && ystart[i]!='2' && ystart[i]!='3' && ystart[i]!='4' && ystart[i]!='5' && ystart[i]!='6' && ystart[i]!='7' && ystart[i]!='8' && ystart[i]!='9' && ystart[i]!='.' && ystart[i]!='E' && ystart[i]!='-' && ystart[i]!='+') {
      printf("A data input ERROR has occured in the X Start tracker position; please go back.\n");
      exit(2);
   }
 }
 for (i=0;i<strlen(expo);i++){
   if (expo[i]!='0' && expo[i]!='1' && expo[i]!='2' && expo[i]!='3' && expo[i]!='4' && expo[i]!='5' && expo[i]!='6' && expo[i]!='7' && expo[i]!='8' && expo[i]!='9' && expo[i]!='.') {
      printf("A data input ERROR has occured in the exposure time; please go back.\n");
      exit(2);
   }
 }
 for (i=0;i<strlen(DEC);i++){
   if(DEC[i]=='+') DEC[i]=' ';
   if(DEC[i]==':') DEC[i]=' ';
   if((DEC[i]=='%' && DEC[i+1]=='3' && DEC[i+2]=='A') || (DEC[i]=='%' && DEC[i+1]=='2' && DEC[i+2]=='B')) {
     DEC[i]=' ';
     if(strlen(DEC)<i+4) { DEC[i+1]=' '; DEC[i+2]=' '; }
       else {
       for(k=i+1;k<strlen(DEC)-2;k++) DEC[k]=DEC[k+2];
       DEC[k]='\0';
     }
   }

   if (DEC[i]!='0' && DEC[i]!='1' && DEC[i]!='2' && DEC[i]!='3' && DEC[i]!='4' && DEC[i]!='5' && DEC[i]!='6' && DEC[i]!='7' && DEC[i]!='8' && DEC[i]!='9' && DEC[i]!='-' && DEC[i]!='.'  && DEC[i]!=' ' && DEC[i]!='+' && DEC[i]!=':') {
      printf("A data input ERROR has occured in the DEC; please go back.\n");
      exit(2);
   }
 }

 sscanf(xstart,"%lf%*c",&xstrt);
 sscanf(ystart,"%lf%*c",&ystrt);
 sscanf(actualaz,"%lf%*c",&obsaz);

 DECdeg = 0;
 DECmin = 0;
 DECsec = 0;

 sscanf(DEC,"%lf %lf %lf%*c",&DECdeg,&DECmin,&DECsec);

 if (DECmin > 60 || DECmin < 0) {
   printf("A data input ERROR has occured in the DEC; please go back.\n");
   exit(2);
 }
 if (DECsec > 60 || DECsec < 0) {
   printf("A data input ERROR has occured in the DEC; please go back.\n");
   exit(2);
 }

 if (DECdeg > 0) DEC2 = DECdeg + DECmin / 60 + DECsec / 3600;
 else DEC2 = DECdeg - DECmin / 60 - DECsec / 3600;

        ha1 = ha2 = -999;

        DEC2 = DEC2/R2D;
        fs = FS;

        az = 1/P*(sin(DEC2)-Q);
        if (fabs(az)>1.0) {
                if (DEC2 > 60.0 /R2D ) az = 0.0;
                else az = 180.0/R2D;
                tde = asin(P * cos(az) + Q);
        } else {
                az = acos(az);
                tde = DEC2;
        }

        c10 = sin(DEC2 - tde);
        c20 = sin(tde);
        c30 = cos(DEC2);
        c40 = c20 * c30;
        c50 = c10 + c40;


        hc = asin ((cos(telangle/R2D)*sin(az))/-cos(tde));
        pa = acos(cos(hc)*cos(az)+sin(hc)*sin(az)*sin(latitude/R2D));
        hc *= R2D / 15;

        inwindow = 0;

        for( ha = -5; ha < 5; ha+=1./60.) {
                h = ha * 15 / R2D ;
                r = -0.1 * fs * cos(DEC2) * sin(h);
                d = fs * (c50 - c20 * cos(DEC2) * cos(h));
                x = r * cos(pa) - d * sin(pa);
                w = r * sin(pa) + d * cos(pa);

                z = fs - (sqrt(fs*fs-x*x-w*w));
                beta = acos(sin(tde) * sin(DEC2) + cos(tde)*cos(DEC2)*cos(h));
                theta = atan(w/(fs*cos(beta)));
                phi = atan(x/(fs*cos(beta)));

                if (inwindow == 0) {
                        if ((((x*x+w*w) < (DXMAX*DYMAX)) &&
                                (x < DXMAX) &&
                                (x > DXMIN) &&
                                (w < DYMAX) &&
                                (w > DYMIN) &&
                                (beta < FBMAX) &&
                                (beta > FBMIN) &&
                                (theta < FTMAX) &&
                                (theta > FTMIN) &&
                                (phi < FPMAX) &&
                                (phi > FPMIN))) {
                        ha1 = ha + hc;
                        inwindow = 1;
                                 }
                 } else if (inwindow == 1){
                         if (!(((x*x+w*w) < (DXMAX*DYMAX)) &&
                                (x < DXMAX) &&
                                (x > DXMIN) &&
                                (w < DYMAX) &&
                                (w > DYMIN) &&
                                (beta < FBMAX) &&
                                (beta > FBMIN) &&
                                (theta < FTMAX) &&
                                (theta > FTMIN) &&
                                (phi < FPMAX) &&
                                (phi > FPMIN))) {
                         ha2 = ha - 1./60. + hc;
                         inwindow = 2;
                         ha = 60;
                         }
                 }
        }
        ha_start=ha1;
        ha_end=ha2;
        az *= R2D;
        if (az == 180.0 || az == 0.0) tracks = 1; else tracks = 2;
        DEC2 = DEC2*R2D;
        if (DEC2 < -12 || DEC2 > 74 ) tracks = 0;
        length = (ha2 - ha1);
        tmidpt = (ha2 + ha1)/2.0;
        hr = length;
        min = (length * 60) - (hr * 60);
        sscanf(expo,"%lf%*c",&exp);
        ttime =  (exp) / 3600.0;

        printf("\n");
        printf("Your target input DEC is %s \n",DEC);
        if (min < 10) printf("This track at Az = %5.1f could have been %d:0%1i (h:mm) long. \n",obsaz,hr,min);
        else printf("This track at Az = %5.1f could have been %d:%2i (h:mm) long. \n",obsaz,hr,min);
        if (obsaz > 58 && obsaz < 78) printf("WARNING: THE CCAS TOWER will obscure some of the mirror AND IS NOT INCLUDED IN THE FILLING FACTOR \n");

        printf("You took a %6.0f second exposure \n",exp);
        printf("\n");

        bmin = 0.0;
        if (tracks == 1 && DEC2 > 0 ) bmin = DEC2 - 65.67;
        else if (tracks == 1) bmin = -(4.3319 + DEC2);
        
        secperdeg = length * 3600 / ((deltatelangle - bmin) * 2.0);

        tstart = - (exp/2.0);
        tend = exp + tstart;

        printf("A perfect exposure of this length would have started %3.0f minutes into the track.\n",( length*30)+tstart/60.0);

     
      fillsum = 0; 
      for(k=tstart;k<=tend;k++) {
         beta = abs(k) / secperdeg  + bmin;
         if (beta <= 1.0 ) fill = 1.0;
         else fill = 1.07 - 0.07 * beta;
         fillsum = fillsum + fill;
/*
         printf ("Integration for %d second at %10.1f beta with filling %10.2f \n",k,beta,fill);
*/
      }
        printf("For this idealized exposure the Average filling factor is %5.2f.\n",fillsum/(exp+1));
        printf("\n");

        betastart = deltatelangle*sqrt(xstrt*xstrt+ystrt*ystrt)/1780;
        if (obsaz > 0 && obsaz < 180 && ystrt > 0 ) betastart = - betastart;
        if (obsaz > 180 && obsaz < 360 && ystrt < 0 ) betastart = - betastart;
        if (az == 180 ) {
           betastart = deltatelangle*sqrt(xstrt*xstrt+(ystrt-bmin/deltatelangle*1780)*(ystrt-bmin/deltatelangle*1780))/1780;
           if (xstrt < 0) betastart = - betastart;
        }
        if (az == 0 ) {
           betastart = deltatelangle*sqrt(xstrt*xstrt+(ystrt-bmin/deltatelangle*1780)*(ystrt-bmin/deltatelangle*1780))/1780;
           if (xstrt > 0 ) betastart = - betastart;
        }
        tstart = betastart * secperdeg;
        tend = tstart + exp;
        fillsum = 0; 
        for(k=tstart;k<=tend;k++) {
           beta = abs(k) / secperdeg  + bmin;
           if (beta <= 1.0 ) fill = 1.0;
           else fill = 1.07 - 0.07 * beta;
           fillsum = fillsum + fill;
        }
        printf("For the actual exposure the Average filling factor is %5.2f.\n",fillsum/(exp+1));

        printf("\n");


        printf("The filling factor ranges from 1.0 (center) to 0.48 (at the edge of the track).\n");
        printf("The maximium filling factor corresponds to a 9.2m aperture with a 3.7m central obscuration:\n");
        printf("     55 square meters of collecting area or a 8.4 meter unobscured telescope \n");

        printf("\n");
        

}