// 
// smeg - A satellite modelling and prediction tool
// Copyright (C) 1999  Tom Rothamel
// 
// 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; see the file COPYING.  If not, write to
// the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
// Boston, MA 02111-1307, USA.

#include <sts2.h>

Locatable *theSun = NULL;

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

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

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

	Sun() {
		init = 0;
	}
};


void Sun::calcpos(double t) {
	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;

	x[0] = px[0] + pdx[0] * t;
	x[1] = px[1] + pdx[1] * t;
	x[2] = px[2] + pdx[2] * t;
}
	
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);
}

void CreateSun() {
	theSun = new Sun(); /* Let there be light! */
}

extern double lat;
extern double lon;
extern double alt;

void ShowSun() {
	Site *site;

	site = new Site(lat, lon, alt);
	
	fprintf(out, "The sun's elevation is currently: %f\n",
		Elevation(theSun, site, sgp_now()));

	delete site;
}

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)

int isLit(Locatable *s, double t) {
	double ox[3];
	double nsx[3];
	double re;
	double rs;
	double q;
	double qe;
	double qs;

	int i;

	theSun->calcpos(t);
	s->calcpos(t);
	
	for (i = 0; i < 3; i++) {
		ox[i] = theSun->x[i] - s->x[i];
		nsx[i] = s->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);

	if (qe > qs && q < qe - qs) {
		return 0;
	} else {
		return 1;
	}
}
	
double SiteSatSun(Locatable *s, Site *si, double t) {
	double sat[3];
	double sun[3];
	int i;

	theSun->calcpos(t);
	s->calcpos(t);
	si->calcpos(t);
	
	for (i = 0; i < 3; i++) {
		sat[i] = s->x[i] - si->x[i];
		sun[i] = theSun->x[i] - si->x[i];
	}

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