Subj:	od
Date:	8/29/00 7:40:00 AM Eastern Daylight Time
From:	Fasd
To:	Fasd

// A minimal Win32 console app

#include <stdio.h>
#include <math.h>
#define pi 3.14159265358979
#define gauss .01720209895
#define obliquity 0.40909280228307

double crossproduct(double x1,double y2,double x2,double y1);
double interpolate(double x1,double x2,double timeinit,double timeobs,double timefinal);
double gaussdays(double julian,double julian2);
double Dectodeg(double deg,double min,double sec);
double RAtodeg(double hr,double min,double sec);
double dotproduct(double x1,double y1,double z1,double x2,double y2,double z2);
double magnitude(double x,double y,double z);
double quadambiguity(double sinambiguously,double cosgayduo);
double radtodeg(double convert);
double degtorad(double convert);
double foftau(double tau,double astdotastdotvector, double rinit);
double goftau(double tau, double rinit);

main() {
/*RA[3][3],/*Dec[3][3]*/
double rhohat[3][3],sunvector[3][3],RAdeci[3],Decdeci[3];
double RArad[3],Decrad[3],tau[3],sunvectorjulian[3],rhohatdotdot[3];
double rhohatdot[3],magsunvector,inclinationec,nodeec,sinperihelionec;
double cosperihelionec,perihelionec,magast2,rhocrossrhodot[3],rhomag;
double deltaast,magrhodot,placeholder[4],meananomaly,eccentricanomaly;
double perihelion,eccentricity,semimajor,magastvector,juliansept,rdotdotrdot;
double rcrossrdot[3],inclination,magrcrossrdot,node,sinnode,cosnode,cosbigU;
double sinbigU,cosnu,sinnu,rdotrdot,astdotvector[3],astvector[3];
double sunvectordot[3],sundotrhohat,nu,bigU,rinitmag,rdotmag, Rmaginit;
double Rmagfinal,goftauforrow2,goftauforrow1,rowvector1[3],rowvector2[3];
double rowmaginit,rowhatfinal[3],rowmagfinal,rowhatinit[3],foftauforrow1;
double foftauforrow2,rvector1[3],rvector2[3],sunvector2[6][3];
double rhocrossrhodotdot[3],astdotastdotvector;
/*juliandate2[3][2],/*UT[3][3]*/

//defining vector components of sunvector from JPL
//sunvector[0]=start hour info given by JPL
//sunvector[1]=info for obs time
//sunvector[2]=end hour info given by JPL
//values surrounding obstimes in order to interpolate sunvector values
/*sunvector2[0][0]=-.6700750614;
sunvector2[0][1]=.6991080949;
sunvector2[0][2]=.3030728689;
sunvector2[1][0]=-.6706015874;
sunvector2[1][1]=.6982431708;
sunvector2[1][2]=.3026979135;
sunvector2[0][0]=-.6700750614;
sunvector2[0][1]=.6991080949;
sunvector2[0][2]=.3030728689;
sunvector2[1][0]=-.6706015874;
sunvector2[1][1]=.6982431708;
sunvector2[1][2]=.3026979135;

UT[0][0]=9;
UT[0][1]=9.38333333;
UT[0][2]=10;

//interpolating x,y,z of first sunvector
sunvector[0][0]=interpolate(sunvector2[0][0],sunvector2[1][0],UT[0][0],
UT[0][1],UT[0][2]);

sunvector[0][1]=interpolate(sunvector2[0][1],sunvector2[1][1],UT[0][0],
UT[0][1],UT[0][2]);
//printf("UT[0][2] %lf",UT[0][2]);
sunvector[0][2]=interpolate(sunvector2[0][2],sunvector2[1][2],UT[0][0],
UT[0][1],UT[0][2]);
//printf("UT[0][2] %lf",UT[0][2]);
sunvector2[2][0]=-.7418705912;
sunvector2[2][1]=.6339141718;
sunvector2[2][2]=.2748075722;
sunvector2[3][0]=-.7308115285;
sunvector2[3][1]=.6448459397;
sunvector2[3][2]=.2795474624;

UT[1][0]=8;
UT[1][1]=8.958333333;
UT[1][2]=9;

//interpolating x,y,z of second sunvector
sunvector[1][0]=interpolate(sunvector2[2][0],sunvector2[3][0],UT[1][0],
UT[1][1],UT[1][2]);
sunvector[1][1]=interpolate(sunvector2[2][1],sunvector2[3][1],UT[1][0],
UT[1][1],UT[1][2]);
sunvector[1][2]=interpolate(sunvector2[2][2],sunvector2[3][2],UT[1][0],
UT[1][1],UT[1][2]);
sunvector2[4][0]=-.8455988726;
sunvector2[4][1]=.5100401014;
sunvector2[4][2]=.2210967209;
sunvector2[5][0]=-.8459808523;
sunvector2[5][1]=.5094931263;
sunvector2[5][2]=.2208595764;

UT[2][0]=9;
UT[2][1]=9.1916666667;
UT[2][2]=10;

//interpolating x,y,z of third sunvector
sunvector[2][0]=interpolate(sunvector2[4][0],sunvector2[5][0],UT[2][0],
UT[2][1],UT[2][2]);
sunvector[2][1]=interpolate(sunvector2[4][1],sunvector2[5][1],UT[2][0],
UT[2][1],UT[2][2]);
sunvector[2][2]=interpolate(sunvector2[4][2],sunvector2[5][2],UT[2][0],
UT[2][1],UT[2][2]);

//julian dates from jpl of hours surrounding obs times from
juliandate2[0][0]=2451759.87500;
juliandate2[0][1]=2451759.91667;
juliandate2[1][0]=2451765.83333;
juliandate2[1][1]=2451765.875;
juliandate2[2][0]=2451775.87500;
juliandate2[2][1]=2451775.91667;

//1st observation data
RA[0][0]=22;
RA[0][1]=1;
RA[0][2]=57.485;
Dec[0][0]=-14;
Dec[0][1]=-11;
Dec[0][2]=-14.33;

sunvectorjulian[0]=interpolate(juliandate2[0][0],juliandate2[0][1],UT[0][0],
UT[0][1],UT[0][2]);
//2nd observation data
RA[1][0]=21;
RA[1][1]=57.17;
RA[1][2]=0;
Dec[1][0]=-14;
Dec[1][1]=-25.9;
Dec[1][2]=0;
sunvectorjulian[1]=interpolate(juliandate2[1][0],juliandate2[1][1],UT[1][0],
UT[1][1],UT[1][2]);

//3rd observation data
RA[2][0]=21;
RA[2][1]=48.83;
RA[2][2]=0;
Dec[2][0]=-14;
Dec[2][1]=-49.999;
Dec[2][2]=0;
sunvectorjulian[2]=interpolate(juliandate2[2][0],juliandate2[2][1],UT[2][0],
UT[2][1],UT[2][2]);

//converting to Decimal degrees
RAdeci[0]=RAtodeg(RA[0][0],RA[0][1],RA[0][2]);
Decdeci[0]=Dectodeg(Dec[0][0],Dec[0][1],Dec[0][2]);
RAdeci[1]=RAtodeg(RA[1][0],RA[1][1],RA[1][2]);
Decdeci[1]=Dectodeg(Dec[1][0],Dec[1][1],Dec[1][2]);
RAdeci[2]=RAtodeg(RA[2][0],RA[2][1],RA[2][2]);
Decdeci[2]=Dectodeg(Dec[2][0],Dec[2][1],Dec[2][2]);

*/

//interpolating x,y,z of first sunvector
//sunvector[1][0]=-.8289400239;

//sunvector[1][1]=.5331247804;
//printf("UT[0][2] %lf",UT[0][2]);
//sunvector[1][2]=.2311054717;
//printf("UT[0][2] %lf",UT[0][2]);
//sunvector2[2][0]=-.7418705912;
//sunvector2[2][1]=.6339141718;
//sunvector2[2][2]=.2748075722;
//sunvector2[3][0]=-.7308115285;
//sunvector2[3][1]=.6448459397;
//sunvector2[3][2]=.2795474624;

//UT[1][0]=8;
//UT[1][1]=8.958333333;
//UT[1][2]=9;

//interpolating x,y,z of second sunvector
//sunvector[0][0]=interpolate(sunvector2[2][0],sunvector2[3][0],UT[1][0],
//UT[1][1],UT[1][2]);
//sunvector[0][1]=interpolate(sunvector2[2][1],sunvector2[3][1],UT[1][0],
/*UT[1][1],UT[1][2]);
sunvector[0][2]=interpolate(sunvector2[2][2],sunvector2[3][2],UT[1][0],
UT[1][1],UT[1][2]);
sunvector2[4][0]=-.8455988726;
sunvector2[4][1]=.5100401014;
sunvector2[4][2]=.2210967209;
sunvector2[5][0]=-.8459808523;
sunvector2[5][1]=.5094931263;
sunvector2[5][2]=.2208595764;

UT[2][0]=9;
UT[2][1]=9.1916666667;
UT[2][2]=10;
*/
//interpolating x,y,z of third sunvector
//sunvector[2][0]=interpolate(sunvector2[4][0],sunvector2[5][0],UT[2][0],
//UT[2][1],UT[2][2]);
//sunvector[2][1]=interpolate(sunvector2[4][1],sunvector2[5][1],UT[2][0],
//UT[2][1],UT[2][2]);
//sunvector[2][2]=interpolate(sunvector2[4][2],sunvector2[5][2],UT[2][0],
//UT[2][1],UT[2][2]);

//julian dates from jpl of hours surrounding obs times from
//juliandate2[0][0]=2451765.83333;
//juliandate2[0][1]=2451765.875;
//juliandate2[2][0]=2451775.87500;
//juliandate2[2][1]=2451775.91667;

//1st observation data
/*RA[1][0]=22;
RA[1][1]=1;
RA[1][2]=57.485;
Dec[1][0]=-14;
Dec[1][1]=-11;
Dec[1][2]=-14.33;*/

//sunvectorjulian[1]=2451773.807;
//2nd observation data
//RA[0][0]=21;
//RA[0][1]=57.17;
//RA[0][2]=0;
//Dec[0][0]=-14;
//Dec[0][1]=-25.9;
//Dec[0][2]=0;
//sunvectorjulian[0]=interpolate(juliandate2[0][0],juliandate2[0][1],UT[0][0],
//UT[0][1],UT[0][2]);

//aditi's entering julian dates
sunvectorjulian[0]=2451765.841;
sunvectorjulian[1]=2451773.807;
sunvectorjulian[2]=2451775.883;

//aditi's entering R for middle observation
sunvector2[2][0]=-.8259730721;
sunvector2[2][1]=.5370605235;
sunvector2[2][2]=.2328119313;
sunvector2[3][0]=-.8263756284;
sunvector2[3][1]=.5365265215;
sunvector2[3][2]=.2325803987;


//aditi's entering interpolated R for all observations

sunvector[0][0]=.7423311584;
sunvector[0][1]=.6334511021;
sunvector[0][2]=.2746067884;
sunvector[1][0]=-.8289400239;
sunvector[1][1]=.5331247804;
sunvector[1][2]=.2311054717;
sunvector[2][0]=.8456720854;
sunvector[2][1]=.5099352645;
sunvector[2][2]=.2210512682;

//3rd observation data
//RA[2][0]=21;
//RA[2][1]=48.83;
//RA[2][2]=0;
//Dec[2][0]=-14;
//Dec[2][1]=-49.999;
//Dec[2][2]=0;
//sunvectorjulian[2]=interpolate(juliandate2[2][0],juliandate2[2][1],UT[2][0],
//UT[2][1],UT[2][2]);

//converting to Decimal degrees
//RAdeci[0]=RAtodeg(RA[0][0],RA[0][1],RA[0][2]);
//Decdeci[0]=Dectodeg(Dec[0][0],Dec[0][1],Dec[0][2]);
//RAdeci[1]=329.3367167;
//Decdeci[1]=-14.76393333;
//RAdeci[2]=RAtodeg(RA[2][0],RA[2][1],RA[2][2]);
//Decdeci[2]=Dectodeg(Dec[2][0],Dec[2][1],Dec[2][2]);

//aditi: entering RA in degrees from Noel
RAdeci[0]=329.33671670;
RAdeci[1]=327.6440542;
RAdeci[2]=327.1901292;
Decdeci[0]=-14.42913333;
Decdeci[1]=-14.76393333;
Decdeci[2]=-14.84758611;


//converting to radians
RArad[0]=degtorad(RAdeci[0]);
Decrad[0]=degtorad(Decdeci[0]);
RArad[1]=degtorad(RAdeci[1]);
Decrad[1]=degtorad(Decdeci[1]);
RArad[2]=degtorad(RAdeci[2]);
Decrad[2]=degtorad(Decdeci[2]);

printf("\n RA in radians: %lf %lf %lf\n", RArad[0], RArad[1], RArad[2]);
printf("\n Dec in radians: %lf %lf %lf\n", Decrad[0], Decrad[1], Decrad[2]);

//convert to gaussian days
tau[0]=gaussdays(sunvectorjulian[0],sunvectorjulian[1]);
tau[1]=gaussdays(sunvectorjulian[2],sunvectorjulian[0]);
tau[2]=gaussdays(sunvectorjulian[2],sunvectorjulian[1]);

printf("\n Tau: %lf %lf %lf \n", tau[0], tau[1], tau[2]);

//finding rhohat for the 3 observations
rhohat[0][0]=cos(RArad[0])*cos(Decrad[0]);
rhohat[0][1]=sin(RArad[0])*cos(Decrad[0]);
rhohat[0][2]=sin(Decrad[0]);
rhohat[1][0]=cos(RArad[1])*cos(Decrad[1]);
rhohat[1][1]=sin(RArad[1])*cos(Decrad[1]);
rhohat[1][2]=sin(Decrad[1]);
rhohat[2][0]=cos(RArad[2])*cos(Decrad[2]);
rhohat[2][1]=sin(RArad[2])*cos(Decrad[2]);
rhohat[2][2]=sin(Decrad[2]);

printf("\n rho hat -1: %lf %lf %lf\n", rhohat[0][0], rhohat[0][1], rhohat[0][2]);
printf("\n rho hat 0: %lf %lf %lf\n", rhohat[1][0], rhohat[1][1], rhohat[1][2]);
printf("\n rho hat 1: %lf %lf %lf\n", rhohat[2][0], rhohat[2][1], rhohat[2][2]);

//printf("\n rhohat -1: %lf\n", rhohat
//solving for rhohatdot

rhohatdot[0]=((pow(tau[2],2)*(rhohat[0][0]-rhohat[1][0]))-(pow(tau[0],2)*(rhohat[2][0]-rhohat[1][0])))/(tau[0]*tau[1]*tau[2]);

rhohatdot[1]=((pow(tau[2],2)*(rhohat[0][1]-rhohat[1][1]))-(pow(tau[0],2)*(rhohat[2][1]-rhohat[1][1])))/(tau[0]*tau[1]*tau[2]);

rhohatdot[2]=((pow(tau[2],2)*(rhohat[0][2]-rhohat[1][2]))-(pow(tau[0],2)*(rhohat[2][2]-rhohat[1][2])))/(tau[0]*tau[1]*tau[2]);

printf("\n rho hat dot: %lf %lf %lf\n ", rhohatdot[0], rhohatdot[1], rhohatdot[2]);


//solving for rhohatdotdot
rhohatdotdot[0]=(-2*((tau[2]*(rhohat[0][0]-rhohat[1][0]))-(tau[0]*(rhohat[2][0]-rhohat[1][0]))))/(tau[0]*tau[1]*tau[2]);

rhohatdotdot[1]=(-2*((tau[2]*(rhohat[0][1]-rhohat[1][1]))-(tau[0]*(rhohat[2][1]-rhohat[1][1]))))/(tau[0]*tau[1]*tau[2]);

rhohatdotdot[2]=(-2*((tau[2]*(rhohat[0][2]-rhohat[1][2]))-(tau[0]*(rhohat[2][2]-rhohat[1][2]))))/(tau[0]*tau[1]*tau[2]);

printf("\n rho hat dot dot: %lf %lf %lf\n ", rhohatdotdot[0], rhohatdotdot[1], rhohatdotdot[2]);

magsunvector=magnitude(sunvector[1][0],sunvector[1][1],sunvector[1][2]);

magastvector=2.5;

//finding cross product between rho and rhodot for each component
rhocrossrhodot[0]=crossproduct(rhohat[1][1],rhohatdot[2],rhohat[1][2],rhohatdot[1]);

rhocrossrhodot[1]=crossproduct(rhohat[1][2],rhohatdot[0],rhohat[1][0],rhohatdot[2]);

rhocrossrhodot[2]=crossproduct(rhohat[1][0],rhohatdot[1],rhohat[1][1],rhohatdot[0]);

//(rhohat cross rhohatdot) dot 
placeholder[0]=dotproduct(rhocrossrhodot[0],rhocrossrhodot[1],rhocrossrhodot[2],sunvector[1][0],sunvector[1][1],sunvector[1][2]);
//printf("rhocrossrhodot[0] %lf rhocrossrhodot[1] %lf rhocrossrhodot [2] %lf sunvector[1][0] %lf sunvector[1][1] %lf sunvector[1][2] %lf",rhocrossrhodot[0], rhocrossrhodot[1] , rhocrossrhodot [2], sunvector[1][0] , sunvector[1][1], sunvector[1][2]);
//(rhohat cross rhohatdot) dot rhohatdotdot
placeholder[1]=dotproduct(rhocrossrhodot[0],rhocrossrhodot[1],rhocrossrhodot[2],rhohatdotdot[0],rhohatdotdot[1],rhohatdotdot[2]);

//printf("placeholder[0] %lf placeholder[1] %lf \n",placeholder[0],placeholder[1]);
printf("\nA is: %lf\n", placeholder[0]/placeholder[1]);
printf("\nB is: %lf\n", ((1+(1/328900.5))/pow(magsunvector,3))*(placeholder[0]/placeholder[1]));


sundotrhohat=dotproduct(sunvector[1][0],sunvector[1][1],sunvector[1][2],rhohat[1][0],rhohat[1][1],rhohat[1][2]);
deltaast=1;

while (deltaast>.0000001 || deltaast<-.0000001) {

rhomag=((1/(pow(magastvector,3)))-((1+(1/328900.5))/pow(magsunvector,3)))*(placeholder[0]/placeholder[1]);
//printf("magastvector %lf",magastvector);
magast2=sqrt(pow(rhomag,2)+pow(magsunvector,2)-(2*sundotrhohat*rhomag));
deltaast=magast2-magastvector;
magastvector+=deltaast;

}
printf("\nrhomag %lf\n",rhomag);
printf("\n r mag %lf \n", magastvector);
 

//rhohat cross rhohatdotdot
rhocrossrhodotdot[0]=crossproduct(rhohat[1][1],rhohatdotdot[2],rhohat[1][2],rhohatdotdot[1]);
rhocrossrhodotdot[1]=crossproduct(rhohat[1][2],rhohatdotdot[0],rhohat[1][0],rhohatdotdot[2]);
rhocrossrhodotdot[2]=crossproduct(rhohat[1][0],rhohatdotdot[1],rhohat[1][1],rhohatdotdot[0]);

placeholder[2]=dotproduct(rhocrossrhodotdot[0],rhocrossrhodotdot[1],rhocrossrhodotdot[2],sunvector[1][0],sunvector[1][1],sunvector[1][2]);
magrhodot=(-.5)*((1/pow(magastvector,3))-((1+(1/328900.5))/pow(magsunvector,3)))*(placeholder[2]/placeholder[1]);
placeholder[3]=placeholder[2]/placeholder[1];

printf("\nrho dot: %lf\n", magrhodot);

//CHECK THIS NEXT STEP
sunvectordot[0]=(sunvector2[3][0]-sunvector2[2][0])/(gauss/24);
sunvectordot[1]=(sunvector2[3][1]-sunvector2[2][1])/(gauss/24);
sunvectordot[2]=(sunvector2[3][2]-sunvector2[2][2])/(gauss/24);

printf("\n R dot vector: %lf %lf %lf\n", sunvectordot[0], sunvectordot[1], sunvectordot[2]);




//Check this one too
astvector[0]=(rhomag*rhohat[1][0])-sunvector[1][0];
astvector[1]=(rhomag*rhohat[1][1])-sunvector[1][1];
astvector[2]=(rhomag*rhohat[1][2])-sunvector[1][2];

printf("\n r vector: %lf %lf %lf\n", astvector[0], astvector[1], astvector[2]);


astdotvector[0]=(magrhodot*rhohat[1][0])+(rhomag*rhohatdot[0])-sunvectordot[0];
astdotvector[1]=(magrhodot*rhohat[1][1])+(rhomag*rhohatdot[1])-sunvectordot[1];
astdotvector[2]=(magrhodot*rhohat[1][2])+(rhomag*rhohatdot[2])-sunvectordot[2];

printf("r dot: %lf %lf %lf", astdotvector[0], astdotvector[1], astdotvector[2]); 
//getting magnitude for rnaught
rinitmag=magnitude(astvector[0],astvector[1],astvector[2]);

//getting magnitude for rnaughtdot
rdotmag=magnitude (astdotvector[0],astdotvector[1],astdotvector[2]);

//getting magnitude for Rinitial
Rmaginit=magnitude (sunvector[0][0],sunvector[0][1],sunvector[0][2]);

//getting magnitude for Rfinal
Rmagfinal=magnitude(sunvector[2][0],sunvector[2][1],sunvector[2][2]);
//finding astvector dot astdotvector for foftau
astdotastdotvector=dotproduct(astvector[0],astvector[1],astvector[2],astdotvector[0],astdotvector[1],astdotvector[2]);

//calculate f of tau for rho init
foftauforrow1=foftau(tau[0],astdotastdotvector,rinitmag);

//calculate g of tau for rho init
goftauforrow1=goftau(tau[0],rinitmag);

//find rvector by combining f of tau and g of tau
rvector1[0]=(foftauforrow1*astvector[0])+(goftauforrow1*astdotvector[0]);
rvector1[1]=(foftauforrow1*astvector[1])+(goftauforrow1*astdotvector[1]);
rvector1[2]=(foftauforrow1*astvector[2])+(goftauforrow1*astdotvector[2]);

//find rowvector by adding rvector and R
rowvector1[0]=rvector1[0]+sunvector[0][0];
rowvector1[1]=rvector1[1]+sunvector[0][1];
rowvector1[2]=rvector1[2]+sunvector[0][2];


//calculate f of tau for row final
goftauforrow2=goftau(tau[1],rinitmag);

//calculate g of tau for row final
foftauforrow2=foftau(tau[1],astdotastdotvector,rinitmag);

//find rvector by combining f of tau and g of tau
rvector2[0]=(foftauforrow2*astvector[0])+(goftauforrow2*astdotvector[0]);
rvector2[1]=(foftauforrow2*astvector[1])+(goftauforrow2*astdotvector[1]);
rvector2[2]=(foftauforrow2*astvector[2])+(goftauforrow2*astdotvector[2]);

//find row vector by adding to R
rowvector2[0]=rvector2[0]+sunvector[2][0];
rowvector2[1]=rvector2[1]+sunvector[2][1];
rowvector2[2]=rvector2[2]+sunvector[2][2];


rowmaginit=magnitude (rowvector1[0],rowvector1[1],rowvector1[2]);
rowmagfinal=magnitude (rowvector2[0],rowvector2[1],rowvector2[2]);
//find final 1st observ row hat
rowhatinit[0]=(rhohat[0][0]/rowmaginit);
rowhatinit[1]=(rhohat[0][1]/rowmaginit);
rowhatinit[2]=(rhohat[0][0]/rowmaginit);

//find final 3rd observ row hat
rowhatfinal[0]=(rhohat[2][0]/rowmagfinal);
rowhatfinal[1]=(rhohat[2][1]/rowmagfinal);
rowhatfinal[2]=(rhohat[2][2]/rowmagfinal);


printf("\n rho hat -1: %lf %lf %lf \n", rowhatinit[0], rowhatinit[1], rowhatinit[2]); 
printf("\n rho hat +1: %lf %lf %lf \n", rowhatfinal[0], rowhatfinal[1], rowhatfinal[2]); 


astdotvector[0]*=gauss;
astdotvector[1]*=gauss;
astdotvector[2]*=gauss;



//finding semimajor axis
rdotdotrdot=dotproduct(astdotvector[0],astdotvector[1],astdotvector[2],astdotvector[0],astdotvector[1],astdotvector[2]);
semimajor=1/((2/magastvector)-(rdotdotrdot/pow(gauss,2)));

printf("\na: %lf\n", semimajor);
//finding eccentricity
//to find h
rcrossrdot[0]=crossproduct(astvector[1],astdotvector[2],astvector[2],astdotvector[1]);
rcrossrdot[1]=crossproduct(astvector[2],astdotvector[0],astvector[0],astdotvector[2]);
rcrossrdot[2]=crossproduct(astvector[0],astdotvector[1],astvector[1],astdotvector[0]);
magrcrossrdot=magnitude(rcrossrdot[0],rcrossrdot[1],rcrossrdot[2]);

eccentricity=sqrt(1-(pow(magrcrossrdot,2)/(pow(gauss,2)*semimajor)));

printf("\ne: %lf\n", eccentricity);
//finding inclination
inclination=sqrt(pow(rcrossrdot[0],2)+pow(rcrossrdot[1],2))/rcrossrdot[2];
inclination=atan(inclination);

printf("\n i: %lf\n", inclination);

//finding long of ascending node
sinnode=rcrossrdot[0]/(magrcrossrdot*sin(inclination));
cosnode=(-1*rcrossrdot[1])/(magrcrossrdot*sin(inclination));

node=quadambiguity(sinnode,cosnode);
printf("\n node %lf \n",node);
//finding perihelion
cosbigU=((astvector[0]*cos(node))+(astvector[1]*sin(node)))/magastvector;
sinbigU=astvector[2]/(magastvector*sin(inclination));

bigU=quadambiguity(sinbigU,cosbigU);

cosnu=(((semimajor*(1-pow(eccentricity,2)))/magastvector)-1)/eccentricity;
rdotrdot=dotproduct(astvector[0],astvector[1],astvector[2],astdotvector[0],astdotvector[1],astdotvector[2]);
sinnu=(semimajor*(1-pow(eccentricity,2))*rdotrdot)/(magastvector*magrcrossrdot*eccentricity);
nu=quadambiguity(sinnu,cosnu);
perihelion=bigU-nu;

printf("\nperihelion %lf \n",perihelion);
//finding eccentric anomaly
eccentricanomaly=(1/eccentricity)*(1-(magastvector/semimajor));
eccentricanomaly=2*pi-acos(eccentricanomaly);
printf("E: %lf", eccentricanomaly);

//finding mean anomaly
juliansept=2451800.5;
meananomaly=eccentricanomaly-(eccentricity*sin(eccentricanomaly));
meananomaly=(gauss/sqrt(pow(semimajor,3)))*(juliansept-sunvectorjulian[1])+meananomaly;
meananomaly=radtodeg(meananomaly);
meananomaly=meananomaly;
printf("\n Mo: %lf\n", meananomaly);


//finding inclination in ecliptic
inclinationec=cos(obliquity)*cos(inclination)+(sin(obliquity)*sin(inclination)*cos(node));
inclinationec=acos(inclinationec);
inclinationec=radtodeg(inclinationec);

printf("\n new i: %lf\n", inclinationec);

//finding long of asc node in ecliptic
sinnode=(sin(node)*sin(inclination))/inclinationec;
cosnode=(cos(obliquity)*cos(inclinationec*pi/180.0)-cos(inclination))/(sin(obliquity)*sin(inclinationec*pi/180.0));

nodeec=quadambiguity(sinnode,cosnode);

nodeec=radtodeg(nodeec);

printf("\n long new: %lf\n", nodeec);

//finding perihelion in ecliptic
sinperihelionec=(sin(obliquity)*sin(node))/sin(inclinationec*pi/180.0);
cosperihelionec=(cos(node)*cos(nodeec*pi/180.0))+(sin(node)*sin(nodeec*pi/180.0)*cos(obliquity));

printf("cos W: %lf", cosperihelionec);
printf("sin W: %lf", sinperihelionec);

perihelionec=quadambiguity(sinperihelionec,cosperihelionec);
perihelionec=-perihelionec+perihelion;
if (perihelionec<0) 
perihelionec=2*pi+perihelionec;
perihelionec=radtodeg(perihelionec);
printf("\n new peri: %lf\n", perihelionec);
}

double degtorad(double convert) {
return convert*(pi/180);
}

double radtodeg(double convert) {
return convert*(180/pi);
}

double quadambiguity(double sinambiguously,double cosgayduo) {
double rightangle;
rightangle=1.57079632679;
	if (sinambiguously>0 && cosgayduo>0) {
	return asin (sinambiguously);
	}
	if (sinambiguously>0 && cosgayduo<0) {
	return (asin(sinambiguously)+rightangle);
	}
	if (sinambiguously<0 && cosgayduo>0) {
	return 2*pi-acos(cosgayduo);
	}
	if (sinambiguously<0 && cosgayduo<0) {
	return (acos(cosgayduo)+rightangle);
	}
}

double magnitude(double x,double y,double z) {
return sqrt(pow(x,2)+pow(y,2)+pow(z,2));
}

double dotproduct(double x1,double y1,double z1,double x2,double y2,double z2) {
return (x1*x2)+(y1*y2)+(z1*z2);
}

double RAtodeg(double hr,double min,double sec) {
return ((hr+(min/60)+(sec/3600))/24)*360;
}

double Dectodeg(double deg,double min,double sec) {
return deg+(min/60)+(sec/3600);
}

double gaussdays(double julian,double julian2) {
return gauss*(julian-julian2);
}

double interpolate(double x1,double x2,double timeinit,double timeobs,double timefinal) {
return x1+(((timeobs-timeinit)/(timefinal-timeinit))*(x2-x1));
}

double crossproduct(double x1,double y2,double x2,double y1) {
return (x1*y2)-(x2*y1);
}

//returns the results of f of tau
double foftau(double tau,double astdotastdotvector, double rinit) {
	return (1-(pow(tau,2)/(2*pow(rinit,3)))+(pow(tau,3)*astdotastdotvector)/(2*pow(rinit,5)));
}

//returns the results of goftau
double goftau(double tau, double rinit) {
return tau-(pow(tau,3)/(6*pow(rinit,3)));
}