/* Compute the lifetime of a plate.  Input is ha, dec, and max offset in
 * arcsec; output is hamin and hamax.  Newer version that adjusts scale
 * as well.
*/

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

static double phi = 32.8;
static double PII=3.141593;
static double R = 58./206265.;

/**********************************************************************/
static double sec2 (double ang) {
   double sec = pow(1./cos(ang),2);
   return sec;
   }

/**********************************************************************/
/* Compute zenith angle given an hour angle and declination
 * in deg. */

double zenith(double h, double dec) {
   double rh, rdec, rphi;
   double rpsi, psi;
   double re, rz, z;
   rh = h*PII/180.;
   rdec = dec*PII/180.;
   rphi = phi*PII/180.;

/* zenith angle */
   re = asin(sin(rdec)*sin(rphi) + cos(rh)*cos(rdec)*cos(rphi));
   rz = PII/2. - re;
   z = rz*180./PII;
   return z;
   }

/**********************************************************************/
/* Compute position angle of north, given an hour angle and dec in deg. */

static double posang(double h, double dec) {
   double rh, rdec, rphi;
   double rpsi, psi;
   double re, rz, z;
   rh = h*PII/180.;
   rdec = dec*PII/180.;
   rphi = phi*PII/180.;

/* Position angle */
   rpsi = atan2(sin(rh)*cos(rphi), cos(rdec)*sin(rphi) - 
	cos(rh)*sin(rdec)*cos(rphi));
   psi = rpsi*180./PII;
   if (psi > 180) psi -= 360.;
   if (psi < -180) psi += 360.;
   return psi;
   }

/******************************************************************/
/* Main routine to compute plate lifetime */
/* Input is ha dec maxoff */

int main (int argc, char **argv) {
   double ha0, dec0, maxoff;
   double za0, rza0, sec2za0, pa0;
   double ha;
   double za, rza, sec2za, pa;
   double theta, drmax0, drmax, ddx, ddy, dr;
   double halim[2], haoff[2];
   int j, i, side;
   double tau, x0[8], y0[8], x[8], y[8];

/* Compute mean change in scale, translation, and rotation, and apply at
 * end - these are factors that are removed by the guider
*/
   double rad0, rad, rot0, rot, rdiff, scale, xavg, yavg, rotavg;
   int n;

   if (argc != 4) {
        printf ("Usage: platelife <ha> <dec> <maxoff>\n");
        exit(1);
        }
   sscanf (argv[1], "%lf", &ha0);
   sscanf (argv[2], "%lf", &dec0);
   sscanf (argv[3], "%lf", &maxoff);

/* Get za, pa at design position */
   pa0 = posang(ha0, dec0);
   za0 = zenith(ha0, dec0);
   rza0 = za0*PII/180.;
   sec2za0 = sec2 (rza0);

   for (j=0; j<2; j++) {
	side = j*2-1;
	halim[j] = ha0;
	haoff[j] = 0;
	for (i=0; i<90; i++) {
	   ha = ha0 + side*i;
	   za = zenith(ha, dec0);
	   rza = za*PII/180.;
	   sec2za = sec2(rza);
	   pa = posang(ha, dec0);
	   theta = (pa-pa0)*PII/180.;

/* Loop through a set of points on the edge of the plate */
	   drmax = 0.;
	   n = 0;
	   xavg = 0;
	   yavg = 0;
	   rotavg = 0;
	   scale = 0;
	   for (tau = 0; tau < 2*PII-.01; tau += PII/4.) {
		x0[n] = sin(tau);
		y0[n] = cos(tau);

/*Convert to arcsec */
		x0[n] = x0[n]*327.*16.5;
		y0[n] = y0[n]*327.*16.5;
		rad0 = sqrt(pow(x0[n],2) + pow(y0[n],2));
		rot0 = atan2(x0[n], y0[n]);

		ddx = (y0[n]*cos(theta)*sin(theta) - x0[n]*pow(sin(theta),2))
		   * pow(tan(rza),2) * R;
		ddy = (-y0[n]*(sec2za - sec2za0) +
		   (x0[n]*sin(theta)*cos(theta)
		   + y0[n]*pow(sin(theta),2))*pow(tan(rza),2)) * R;

		x[n] = x0[n] + ddx;
		y[n] = y0[n] + ddy;
		rad = sqrt(pow(x[n],2) + pow(y[n],2));
		rot = atan2(x[n], y[n]);
		rdiff = rot-rot0;
		if (rdiff < -PII) rdiff += 2.*PII;
		if (rdiff > PII) rdiff -= 2.*PII;
		rotavg += rdiff;
		xavg += (x[n]-x0[n]);
		yavg += (y[n]-y0[n]);
		scale += (rad/rad0-1.);
		n++;
		}
	   xavg = xavg/n;
	   yavg = yavg/n;
	   rotavg = rotavg/n;
	   scale = scale/n;

/* Reset drmax */
	   drmax = 0.;
	   n=0;
	   for (tau = 0; tau < 2*PII; tau += PII/4.) {

		rad = sqrt(pow(x[n],2) + pow(y[n],2));
		rot = atan2(x[n], y[n]);

/* Remove mean corrections */
		rad = rad/(1.+scale);
		rot = rot - rotavg;
		x[n] = rad*sin(rot)-xavg;
		y[n] = rad*cos(rot)-yavg;
		ddx = x[n]-x0[n];
		ddy = y[n]-y0[n];

		dr = sqrt(pow(ddx,2) + pow(ddy,2));
		if (dr > drmax) drmax = dr;
		n++;
		}
	   if (drmax > maxoff) break;
	   halim[j] = ha;
	   haoff[j] = i * side;
	   }
	}
   printf ("{east %.0f} {west %.0f}\n", halim[0], halim[1]);

   return 0;
   }