00021 {
00022
00023 using namespace std;
00024
00025
00026 double ra=0;
00027 double altitude = 550.e3 ;
00028 double incl = 28.5;
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
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
00067 double duration = SimDurationTime;
00068
00069 long Nsteps = (long)(duration/delt+0.5) + 1L;
00070
00071 double Tlive=0.;
00072
00073
00074 double elapse = 0.;
00075
00076
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
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
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
00178
00179 #ifdef comparison
00180
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;
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 }
00233
00234 out.close();
00235 return 0;
00236 }
