Orbit.cxx Source File

00001 
00002 #include <string>
00003 #include <cmath>
00004 #include <iostream>
00005 #include <fstream> 
00006 
00007 #include "astroutil.h"
00008 #include "solsystem.h"
00009 
00010 #define comparison
00011 #ifdef comparison // for comparison
00012 #include "astro/JulianDate.h"
00013 #include "astro/EarthOrbit.h"
00014 #include "astro/SolarSystem.h"
00015 #include "astro/EarthCoordinate.h"
00016 #include <cassert>
00017 #endif
00018 #define zenmax  105.F
00019 #define Ao      13000.F // cm^2 maximum area at 0 angle incidence
00020 
00021 main(){
00022     
00023     using namespace std;
00024     //source coordinates
00025     
00026     double ra=0;
00027     double altitude = 550.e3  ; //m
00028     double incl = 28.5; // deg
00029     double e = 0. ;
00030     double rock = 35.; 
00031     double delt = 30.; 
00032     double delthmin = 5.;
00033     double delthmax = 10.;
00034     double ddec=2.0;
00035     long Norb=0;
00036     
00037     const double TliveEff=0.90;
00038     //Time
00039     const double MissionStartTime =GetJD(2005,1,1,0.0);
00040     double JDStart=GetJD(2005.,7,18,0.0);
00041     double StartSimDate=(GetJD(2005.,7,18,0.0)-MissionStartTime)*SecondsPerDay;
00042     double GMST0=GetGMST(GetJD(2005.,7,18,0.0));
00043     double EndSimDate=(GetJD(2005,7,19,0.0)-MissionStartTime)*SecondsPerDay;
00044     double SimDurationTime =(EndSimDate-StartSimDate);
00045     
00046     
00047     double sini = sin(incl*D2R);
00048     double cosi = cos(incl*D2R);
00049     
00050     
00051     double Rearth = R_TERRA  ; 
00052     double alt=(Rearth + altitude)/1000.;
00053     double a = (Rearth + altitude)/Rearth; 
00054     
00055     double T = PI2*pow(a,1.5)/sqrt(5.98e24*6.67e-11)*pow(Rearth,1.5)  ;         
00056     
00057     double u = 3.9860044e14/pow(Rearth,3)  ; 
00058     double n = sqrt(u / pow(a,3));
00059     double n1 = n*(1. + 3.*J2/2.*sqrt(1. - e*e)/(a*a)/pow((1. - e*e),2)*(1.-1.5*sini*sini));
00060     double dwdt = n1*3.*J2/2./(a*a)/pow((1. - e*e),2)*(2. - 2.5*sini*sini);
00061     
00062     double dOmegadt = -n1*3.*J2/2./(a*a)/pow((1. - e*e),2)*cosi;
00063     
00064     double dMdt = n1;
00065     
00066     //double duration = 2.*PI/fabs(dOmegadt); 
00067     double duration = SimDurationTime;
00068     
00069     long Nsteps = (long)(duration/delt+0.5) + 1L;
00070     
00071     double Tlive=0.;
00072     
00073     // Initialize the elapsed time
00074     double elapse = 0.;
00075     
00076     // Convergence tolerance for determining the eccentric anomaly below
00077     int IsSAA = 0; 
00078     
00079     double M0=dMdt*StartSimDate;
00080     range(&M0,6.28);
00081     double Omega0 = dOmegadt*StartSimDate;
00082     double w0 = dwdt*StartSimDate;
00083     double M,Omega,w;
00084     double StartDate, EndDate;
00085     //Nsteps=4;
00086     ofstream out("orbit2.out");
00087     SolSystem Sun, Moon;
00088     Sun.SetLocation(0.,0.,0.);
00089     Moon.SetLocation(0.,0.,0.);
00090     double SunRa;//
00091     double SunDec;//
00092     double MoonRa;//
00093     double MoonDec;//
00094     double longitude;
00095     double latitude;
00096     double LATraZ;
00097     double LATDecZ;
00098     double LATraX;
00099     double LATDecX;
00100     double SCPos[3];
00101     double SCLat,SCLon,SCAlt;
00102     double RaZ,DecZ;
00103 
00104 #ifdef comparison
00105     astro::EarthOrbit orbit; 
00106 #endif
00107 
00108     for (long ti = 0;ti<(Nsteps-1);ti++){
00109         if( ti%100 ==0) cout << ti << endl;
00110         elapse = elapse + delt;
00111         StartDate=StartSimDate+elapse;
00112         JDStart+=(delt/86400.);
00113         GMST0=GetGMST(JDStart);
00114         EndDate=StartDate+delt;
00115         
00116         Sun.SetObj(SUN);
00117         Sun.CalculatePos(JDStart);
00118         Moon.SetObj(MOON);
00119         Moon.CalculatePos(JDStart);
00120 
00121         SunRa=Sun.GRa;
00122         SunDec=Sun.GDec;
00123         MoonRa=Moon.GRa;
00124         MoonDec=Moon.DDec;
00125         
00126         M=M0+dMdt*elapse;
00127         range(&M,6.28);
00128         Omega = Omega0+dOmegadt*elapse;
00129         range(&Omega,6.28);
00130         w = w0+dwdt*elapse;
00131         double Enew;
00132         
00133         Enew=Kepler(M,e);
00134         
00135         double pp[3],location[3],p[3][3],ppr[3],pointdir[3],zenitdir[3];
00136         double a1[3] = {a*(cos(Enew) - e),a*sqrt(1.- e*e) * sin(Enew),0};
00137         calc_unit_vector(a1, pp);
00138         double cosOmega = cos(Omega)  ; 
00139         double sinOmega = sin(Omega);
00140         double cosw = cos(w) ;         
00141         double sinw = sin(w);
00142         
00143         p[0][0] = cosw*cosOmega - sinw*sinOmega*cosi;
00144         p[0][1] = cosw*sinOmega + sinw*cosOmega*cosi;
00145         p[0][2] = sinw*sini;
00146         p[1][0] = -sinw*cosOmega - cosw*sinOmega*cosi;
00147         p[1][1] = -sinw*sinOmega + cosw*cosOmega*cosi;
00148         p[1][2] = cosw*sini;
00149         p[2][0] = sinOmega*sini;
00150         p[2][1] = -cosOmega*sini;
00151         p[2][2] = cosi;
00152         
00153         vector_matrix_multiply(pp,p,location);
00154         memcpy(zenitdir,location,3*sizeof(double));
00155         SCPos[0]= zenitdir[0]*alt;
00156         SCPos[1]= zenitdir[1]*alt;
00157         SCPos[2]= zenitdir[2]*alt;
00158         
00159         longitude=atan2(location[1],location[0])*R2D;
00160         RaZ=longitude;
00161         if(RaZ<0.)RaZ+=360;
00162         //longitude = (longitude +(0.25068447/60.*elapse*SidSolar) +(GMST0+elapse*SidSolar*15./3600.); //verificare 1.00273790935
00163         longitude=(GMST0-RaZ);
00164         if(longitude>180.)longitude-=360.;
00165         if(longitude<-180.)longitude+=360;
00166         
00167         latitude=asin(location[2]);
00168         DecZ=latitude*R2D;
00169         latitude=(atan(tan(latitude))/((1.-EarthFlat)*(1.-EarthFlat)) );
00170         
00171         SCAlt=sqrt(SCPos[0]*SCPos[0]+SCPos[1]*SCPos[1])/cos(latitude)-Rearth/(1000.*sqrt(1.-0.00669454*0.00669454*sin(latitude)*sin(latitude)));
00172         SCLat=latitude*R2D;
00173         Tlive=delt*TliveEff;
00174         IsSAA=InsideSAA(longitude,latitude);
00175         SCLon=longitude;
00176         int temp = (int)(elapse/T); 
00177         //range(&longitude,360.);
00178    
00179 #ifdef comparison
00180         // temporary to compare with copy
00181         Hep3Vector pos = orbit.position(astro::JulianDate(JDStart));
00182         Hep3Vector posdiff(pos.x()-SCPos[0], pos.y()-SCPos[1], pos.z()-SCPos[2]);
00183         double absdiff=posdiff.mag();
00184 
00185         astro::EarthCoordinate epos(pos, JDStart);
00186         double lat=epos.latitude(), lon=epos.longitude(), alt=epos.altitude();
00187 
00188         astro::SkyDir sun = astro::SolarSystem(astro::SolarSystem::Sun, JDStart );
00189         astro::SkyDir moon = astro::SolarSystem(astro::SolarSystem::Moon, JDStart );
00190 
00191         double sun_ra = sun.ra(), sun_dec= sun.dec();
00192         double moon_ra = moon.ra(), moon_dec = moon.dec();
00193 
00194         bool test =  absdiff< 0.050 && fabs(lat-SCLat) < 0.2 && fabs(lon-SCLon)<0.2 ;
00195         if( !test){
00196             cout << "failed: elapse=" << elapse << endl;
00197         }
00198 #endif
00199         GetRockMat(pp,rock,temp,ppr);
00200         vector_matrix_multiply(ppr,p,pointdir);
00201         
00202         LATraZ=atan2(pointdir[1],pointdir[0])*R2D;
00203         LATraZ = LATraZ; 
00204         range(&LATraZ,360.);
00205         LATDecZ=asin(pointdir[2])*R2D;
00206         
00207         LATDecX=LATDecZ; //+0.25068447/60.*elapse* +GMST0;
00208         LATraX=LATraZ-90;
00209         range(&LATraX,360.);
00210         out.flags(ios::fixed);
00211         
00212         out<<StartDate <<'\t';
00213         out<<EndDate <<'\t';
00214         out<< SCPos[0]<<'\t';
00215         out<< SCPos[1]<<'\t';
00216         out<< SCPos[2]<<'\t';
00217         out<<LATraZ<<'\t';
00218         out<<LATDecZ<<'\t';
00219         out<<LATraX<<'\t';
00220         out<<LATDecX<<'\t';
00221         out<<RaZ<<"\t";
00222         out<<DecZ <<'\t';
00223         out<<"1"<<'\t';
00224         out<<Tlive <<'\t';
00225         out<<IsSAA<<'\t';
00226         out<<SCLon <<'\t';
00227         out<<SCLat <<'\t';
00228         out<<SCAlt <<'\t';
00229         out<< SunRa <<"\t"<<SunDec  <<'\t';
00230         out<<MoonRa <<"\t"<<MoonDec   <<endl;
00231         
00232                                 }//time loop
00233                                 
00234                                 out.close();
00235                                 return 0;
00236 }
00237 
00238 
00239 
00240 
00241 //Contents Units
00242 ;
00243 /*out<<"1 Start (Mission Elapsed Time) s: "<<StartDate <<endl;cout<<"2 End (Mission Elapsed Time) s :"<<EndDate <<endl;
00244 out<<"3 position of S/C at start of interval (x,inertial coordinates) km "<<  zenitdir[0]*alt/1000. <<endl;
00245 out<<"3 position of S/C at start of interval (y,inertial coordinates) km "<<  zenitdir[1]*alt/1000.<<endl;
00246 out<<"3 position of S/C at start of interval (z,inertial coordinates) km "<<  zenitdir[2]*alt/1000.<<endl;
00247 //cout<<"4 viewing direction at start (LAT +z axis), 2 angles deg"<< <<endl;
00248 //cout<<"5 orientation at start (LAT +x axis), 2 angles deg"<< <<endl;
00249 out<<"6 zenith direction at start, 2 angles deg "<<" "<<longitude<<" "<<latitude <<endl;
00250 out<<"7 LAT operation mode dimensionless: "<<"1"<<endl;
00251 out<<"8 livetime s:"<<Tlive <<endl;
00252 out<<"9 SAA flag dimensionless :"<<saa_count<<endl;
00253 out<<"10 S/C longitude deg :"<<longitude <<endl;
00254 out<<"11 S/C latitude deg :"<<latitude <<endl;
00255 out<<"12 S/C altitude km :"<<alt/1000. <<endl;
00256 out<<"13 direction of the sun, 2 angles Deg :"<< SunRa <<" "<<SunDec  <<endl;
00257 out<<"14 direction of the moon, 2 angles deg: "<<MoonRa <<" "<<MoonDec   <<endl;*/
00258