// sattool - visual satellite tracking and prediction tool.
// Copyright 2000 Tom Rothamel <tom-idbg@onegeek.org>
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  

#include "sattool.h"
#include "vector.h"

class Sun : public Object {
public:	
	int init;
	double ct; /* Start of 2 day interpolation period. */

	double px[3];
	double pdx[3];

	void calcrawpos(double *, double);
	Position *calcpos(double, Position *);

	Sun() {
		init = 0;
		lock(this);
	}
};

Sun *theSun = NULL;

Position *Sun::calcpos(double t, Position* p) {
	if (!p) p = new Position;
	t -= ct;
	
	if (!init || t < 0 || t > 2 ) {
		init = 1;

		ct += t - 1;
		
		calcrawpos(px, ct);
		calcrawpos(pdx, ct + 2);

		pdx[0] -= px[0];
		pdx[1] -= px[1];
		pdx[2] -= px[2];

		t = 1;
	}
	
	t /= 2;

	p->x[0] = px[0] + pdx[0] * t;
	p->x[1] = px[1] + pdx[1] * t;
	p->x[2] = px[2] + pdx[2] * t;

	return p;
}
	
double rev( double x ) {
        return x - floor(x/360.0)*360.0;
}

#define sind(x) sin((x) * DEGRAD)
#define cosd(x) cos((x) * DEGRAD)
#define atan2d(x, y) (atan2(x, y) * RADEG)

/* This code sucks. It needs to be rewritten so that it doesn't keep
   converting to and from degrees. I'm also not sure anymore how
   accurate it is.
*/
void Sun::calcrawpos(double *lpx, double t) {
	double w, a, e, M, oblecl, L, E, qx, qy, qz, r, lon, v, d;
	
	d = t;
	w = 282.9404 + 4.70935E-5 * d;
	a = 149598073; /* 1.0 AU; */
	e = 0.016709 - 1.151E-9 * d;
	M = 356.0470 + 0.9856002585 * d;
	M = rev(M);
	
        E = M + (180/M_PI) * e * sind(M) * (1 + e * cosd(M));
        
        qx = a * (cosd(E) - e);
        qy = a * sind(E) * sqrt(1 - e * e);

        r = sqrt(qx*qx + qy*qy);
        v = atan2d(qy, qx);

        lon = rev(v + w);

        qx = r * cosd(lon);
        qy = r * sind(lon);

	lpx[0] = qx;
	lpx[1] = qy * cosd(23.4406);
	lpx[2] = qy * sind(23.4406);
}

double arcsin(double arg) {
	return atan(arg/sqrt(1-(arg*arg)));
}

double arccos(double arg) {
	return atan2(sqrt(1-arg*arg), arg);
}

#define SR 696000.0 /* Solar Radius */

#define show(x) fprintf(stderr, "%s = %f\n", #x, x)

// This determines if a position is lit or not. This only takes the earth
// and the sun into account. (So something in the moon's shadow would
// still be considered lit.)
int isLit(Position *p) {
	Position *sunpos;
	double t = p->t;
	double ox[3];
	double nsx[3];
	double re;
	double rs;
	double q;
	double qe;
	double qs;

	int i;

	if (!theSun) theSun = new Sun;
      	sunpos = theSun->calcpos(t, NULL);

	for (i = 0; i < 3; i++) {
		ox[i] = sunpos->x[i] - p->x[i];
		nsx[i] = p->x[i] * -1;
	}
	
	re = len(nsx); /* Sat -> Earth */
	rs = len(ox); /* Sat -> Sun */

	q = dot(nsx, ox)/(re * rs);
	q = arccos(q);

	qe = arcsin(R/re);
	qs = arcsin(SR/rs);

	delete sunpos;
	
	if (qe > qs && q < qe - qs) {
		return 0;
	} else {
		return 1;
	}
}
	
// This figures out the Site--Satellite--Sun angle. (Used to calculate
// the fraction illuminated.)
double SiteSatSun(Position *s, Position *si) {
	Position *sunpos;
	double t = s->t;
	double sat[3];
	double sun[3];
	int i;

	if (!theSun) theSun = new Sun;
	sunpos = theSun->calcpos(t, NULL);
	
	for (i = 0; i < 3; i++) {
		sat[i] = s->x[i] - si->x[i];
		sun[i] = sunpos->x[i] - si->x[i];
	}

	delete sunpos;
	
	return fabs(arccos(dot(sat, sun)/(len(sat)*len(sun))));
}

// This returns the sun elevation at the given site at the given time.
double SunElevation(Object *site, double t) {
	double el;
	Encounter *e;

	if (!theSun) theSun = new Sun;

	e = new Encounter(theSun, site, t);
	el = e->elevation;
	delete e;

	return el;
}