/*	determines the time cost (in seconds ) of a given move
	Note:

		1)	cost =	move_time+penalties
		2)  the lst is updated for the move, delay and exposure time
*/

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

#define TELESCOPE_SLEW_RATE		1.0		/* in degrees per second */
#define TELESCOPE_SETTLE_TIME	1.0		/* in seconds 			*/
#define TELESCOPE_LATITUDE		38.0	/* in degrees 			*/
#define TELESCOPE_DOME_PENALTY	15.0	/* zing em if they move it */
#define TELESCOPE_DOME_RATE		2.0		/* degees/per second 	*/ 
#define TELESCOPE_DOME_SETTLE	4.0

#define EXTINCTION				.1		/* magnitudes per airmass */
#define READOUT_TIME			10.0	/* CCD readout time 	*/
#define BIG_SECZ				50.0

static double move_time(double lst,TARGET from ,TARGET to,double *penalty);
static double obs_time(double lst,TARGET to, double *penalty);
static double azimuth(double, double, double *);
static double bright_to_time(double bright, double mag);
static void check_coords(TARGET);

static double extinction=EXTINCTION;

extern verbose;

double cost(lst, from, to)
double *lst;
TARGET from, to;
{
	double t_move,t_obs,p_move,p_obs,cost;

	check_coords(from);
	check_coords(to);

	t_move= move_time(*lst,from,to,&p_move);
	t_obs = obs_time(*lst,to,&p_obs) ;
	(*lst)+=(t_move + t_obs + to.delay )*SECONDS_TO_DEGREES;

	cost = t_move + p_move + p_obs;

	if(verbose & 2){
		static done;
		if(!done){
			printf("(ra,dec)->(ra,dec)   cost  tmv   p_move  p_obs\n");
			done=1;
		}
		printf("(%7.2f,%7.2f) -> (%7.2f,%7.2f) %7.1f %7.1f %7.1f %7.1f\n",
			from.ra,from.dec,
			to.ra,to.dec,
			cost,t_move,p_move,p_obs);
	}
	return cost;
}

void modify_extinction(in)
double in;
{
	extinction=in;
}

/*	computes the time to move the telescope between
	objects.  The time is the maximum of the telescope
	motion and dome times.  The dome will not move
	until the azimuth is greater than the tolerance function.
	Provision for increasing the cost of a dome move is included
*/

static double move_time(lst,from,to,penalty)
double lst,*penalty;
TARGET from,to;
{
	double dra,ddec,dist,tel_time,dome_time,dome_tol,dome_move;
	double from_az,to_az,secz;

	dra=fabs(fmod2(from.ra-to.ra,360.));		/* calculate the telescope move time */
	ddec=fabs(from.dec-to.dec);
	dist=(dra>ddec) ? dra : ddec ;
	tel_time=dist/TELESCOPE_SLEW_RATE+TELESCOPE_SETTLE_TIME;

	from_az=azimuth(lst-from.ra,from.dec,&secz);
	to_az=azimuth(lst-to.ra,to.dec,&secz);
	dome_move=fabs(fmod2(from_az-to_az,360.));
	dome_tol=5.0/sqrt(1.001-1./pow(secz,2.0));		/* whether to move slit */
	if(dome_move>dome_tol){
		dome_time=dome_move/TELESCOPE_DOME_RATE+TELESCOPE_DOME_SETTLE;
		*penalty=TELESCOPE_DOME_PENALTY;
	}else{
		dome_time=0.;
		*penalty=0.;
	}

	return (tel_time>dome_time) ? tel_time : dome_time;
}

/*	returns the extra time required to observe a target 
	due to the being off the meridian
		1. due to extinction
	MOON has not been included
	TWILITE
*/

static double obs_time(lst,in,penalty)
double lst,*penalty;
TARGET in;
{
	double secz,best_secz,t_ext,t_moon,t_sky,rval;

	best_secz=1.0/cos(PF*(TELESCOPE_LATITUDE-in.dec));
	azimuth(lst-in.ra,in.dec,&secz);
	if(secz>=1.0){
		t_ext=in.exposure*pow(10.0,2.0*extinction*(secz-best_secz)/2.5);
		t_moon=bright_to_time(moon_light(lst,in.ra,in.dec),in.mag);
		t_sky=bright_to_time(sky_light(lst,in.ra,in.dec),in.mag);
		*penalty= t_ext + t_moon + t_sky ;
	 	rval=in.exposure + READOUT_TIME ;
	}else{
		*penalty=BIG_TIME;
		rval=0.;				/* object will not be done */
	}		
	return rval;
}

/*	calculates extra time required to reach magnitude limit
	due to bright background
*/

static double bright_to_time(bright,mag)
double bright,mag;
{

	if(mag<bright)				
		return 0. ;
	else
		return BIG_TIME;
}	

static double azimuth(ha,dec,psecz)
double ha,dec,*psecz;
{
	static double sl,cl,latitude;
	double az,sh,sd,ch,cd,y,temp;

	if(latitude==0.){				/* calculate for once and for all */
		latitude=TELESCOPE_LATITUDE;
		sl=sin(latitude*PF);
		cl=cos(latitude*PF);
	}

	sh=sin(ha*PF);	sd=sin(dec*PF);	
	ch=cos(ha*PF);	cd=cos(dec*PF);	

	temp=sl*sd+cl*cd*ch;
	if(temp>(1.0/BIG_SECZ))
		*psecz=1.0/temp;
	else
		*psecz=BIG_SECZ;

	y=ch*sl-sd*cl/cd;				/*calc telescope azimuth */

	if((fabs(sh)+fabs(y))<SMALL)	/* check for zenith condition */
		az=180.;
	else
		az=atan2(sh,y)/PF+180.;
	return(az);
}

static void check_coords(in)
TARGET in;
{
	if( in.dec<-90. || in.dec>90.)
		error(1,"invalid dec");

	if( in.ra<0. || in.ra>360.)
		error(1,"invalid  ra");
}