// 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 <malloc.h>
#include <string.h>
#include <unistd.h>

// An EarthSite is a site that's fixed to the surface of the earth.
EarthSite::EarthSite(char *n, double lat, double lon, double altin) {
	int nl;

	storelat = lat;
	storelon = lon;
	
	sinlat = sin(lat * DEGRAD);
	coslat = cos(lat * DEGRAD);

	lonr = lon * DEGRAD;
	alt = altin;

	name = strdup(n);
}

void EarthSite::write(FILE *f, int d) {
	pad(f, d);
	fprintf(f, "(earthsite \"%s\" %f %f %f)\n", name, storelat,
		storelon, alt);
}

EarthSite::~EarthSite() {
	free(name);
}

static const double F = (1.0/298.257222);
static const double sqomf = (1.0-F)*(1.0-F); 
static const double ffff = (F + F - F * F);
static const double ROTRATE = ((2*M_PI/86400.0)*1.002737909350); 

Position *EarthSite::calcpos(double t, Position *p) {
	double c;
	double s;
	double theta;

	if (!p) p = new Position;

	theta = sgp_gmst(t) + lonr;

	p->t = t;
	p->coslat = coslat;
	p->sinlat = sinlat;
	p->costheta = cos(theta);
	p->sintheta = sin(theta);
	
	c = 1.0/sqrt(1 - ffff * sinlat * sinlat);
	s = sqomf * c;

 	p->x[0] = (R*c + alt) * p->coslat * p->costheta;
	p->x[1] = (R*c + alt) * p->coslat * p->sintheta;
	p->x[2] = (R*s + alt) * p->sinlat;

	p->dx[0] = -ROTRATE * p->x[1];
	p->dx[1] = ROTRATE * p->x[0];
	p->dx[2] = 0;

	p->sinlat *= -1;
	p->sintheta *= -1;

	return p;
}