/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*/final double E6A = 1.E-6;final double TOTHRD = .66666667;final double XJ2 = 1.082616E-3;final double XJ4 = -1.65597E-6;final double XJ3 = -.253881E-5;final double XKE = .743669161E-1;final double XKMPER = 6378.135;final float erad = 6378.135;final double XMNPDA = 1440.;final double AE = 1.;final double QO = 120.0;final double SO = 78.0;final double CK2 =.5*XJ2*Math.pow(AE,2);final double CK4 = -.375*XJ4*Math.pow(AE,4);final double QOMS2T = Math.pow(((QO-SO)*AE/XKMPER),4);final double S = AE*(1.+SO/XKMPER);final static double DE2RA = .174532925E-1;final static double PI = 3.14159265;final static double PIO2 = 1.57079633;final static double TWOPI = 6.2831853;final static double X3PIO2 = 4.71238898;
class Camera{float cx,cy,cz,tx,ty,tz;float shotLength;public Camera(float cx,float cy, float cz, float tx, float ty, float tz){this.cx = cx;this.cy = cy;this.cz = cz;this.tx = tx;this.ty = ty;this.tz = tz;calcShotLength();}void calcShotLength(){shotLength = sqrt((cx - tx)*(cx-tx) + (cy - ty)*(cy-ty) + (cz - tz)*(cz-tz));}void set(float x, float y, float z){cx = x;cy = y;cz = z;}void feed(){//left handed coordinate system pleasecamera(cx,cy,cz,tx,ty,tz,0,0,1);scale(-1,-1,-1);}void rotZ(float rads){float tcx = cx;cx = cx*cos(rads) + cy*sin(rads);cy = -tcx*sin(rads) + cy*cos(rads);//rounding errors are not energy preserving//normalizefloat len = sqrt((cx - tx)*(cx-tx) + (cy - ty)*(cy-ty) + (cz - tz)*(cz-tz));cx /= len;cy /= len;cz /= len;cx *= shotLength;cy *= shotLength;cz *= shotLength;}void upZ(float amt){//this is weird :/cz += shotLength*amt;if (cz > shotLength) cz = shotLength;if (cz < -shotLength) cz = - shotLength;float len = sqrt((cx - tx)*(cx-tx) + (cy - ty)*(cy-ty) + (cz - tz)*(cz-tz));cx /= len;cy /= len;cz /= len;cx *= shotLength;cy *= shotLength;cz *= shotLength;}void modShot(float amt){float len = getLength();shotLength *= amt;if (camMode == 0){if (shotLength < erad + 100) shotLength = erad + 100;}else{if (shotLength < 100) shotLength = 100;}if (shotLength > erad + 120000) shotLength = erad + 120000;cx /= len;cy /= len;cz /= len;cx *= shotLength;cy *= shotLength;cz *= shotLength;}float getLength(){return sqrt((cx - tx)*(cx-tx) + (cy - ty)*(cy-ty) + (cz - tz)*(cz-tz));}float getHeight(){return sqrt((cx)*(cx) + (cy )*(cy) + (cz)*(cz));}}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*/public class DEEP {double DAY,PREEP,XNODCE,ATIME,DELT,SAVTSN,STEP2,STEPN,STEPP;double ZNS;double C1SS;double ZES;double ZNL;double C1L;double ZEL;double ZCOSIS;double ZSINIS;double ZSINGS;double ZCOSGS;double ZCOSHS;double ZSINHS;double Q22;double Q31;double Q33;double G22;double G32;double G44;double G52;double G54;double ROOT22;double ROOT32;double ROOT44;double ROOT52;double ROOT54;double THDT;{ZNS = 1.19459E-5f;C1SS = 2.9864797E-6f;ZES = .01675f;ZNL = 1.5835218E-4f;C1L = 4.7968065E-7f;ZEL = .05490f;ZCOSIS = .91744867f;ZSINIS = .39785416f;ZSINGS = -.98088458f;ZCOSGS = .1945905f;ZCOSHS = 1.0f;ZSINHS = 0.0f;Q22 = 1.7891679E-6f;Q31 = 2.1460748E-6f;Q33 = 2.2123015E-7f;G22 = 5.7686396f;G32 = 0.95240898f;G44 = 1.8014998f;G52 = 1.0508330f;G54 = 4.4108898f;ROOT22 = 1.7891679E-6f;ROOT32 = 3.7393792E-7f;ROOT44= 7.3636953E-9f;ROOT52 = 1.1428639E-7f;ROOT54 = 2.1765803E-9f;THDT = 4.3752691E-3f;}double THGR;double EQ;double XNQ;double AQNV;double XQNCL;double XMAO;double XPIDOT;double SINQ;double COSQ;double OMEGAQ;double STEM;double CTEM;double ZCOSIL;double ZSINIL;double ZSINHL;double ZCOSHL;double C;double GAM;double ZMOL;double ZX;double ZY;double ZCOSGL;double ZSINGL;double ZMOS;int LS;double ZCOSG;double ZSING;double ZCOSI;double ZSINI;double ZCOSH;double ZSINH;double CC;double ZN;double ZE;double ZMO;double XNOI;double A1;double A3;double A7;double A8;double A9;double A10;double A2;double A4;double A5;double A6;double X1;double X2;double X3;double X4;double X5;double X6;double X7;double X8;double Z31;double Z32;double Z33;double Z1;double Z2;double Z3;double Z11;double Z12;double Z13;double Z21;double Z22;double Z23;double S1;double S2;double S3;double S4;double S5;double S6;double S7;double SE;double SI;double SL;double SGH;double SH;double EE2;double E3;double XI2;double XI3;double XL2;double XL3;double XL4;double XGH2;double XGH3;double XGH4;double XH2;double XH3;double SSE;double SSI;double SSL;double SSH;double SSG;double SE2;double SI2;double SL2;double SGH2;double SH2;double SE3;double SI3;double SL3;double SGH3;double SH3;double SL4;double SGH4;int IRESFL;int ISYNFL;double EOC;double G201;double G211;double G310;double G322;double G410;double G422;double G520;double G533;double G521;double G532;double SINI2;double F220;double F221;double F321;double F322;double F441;double F442;double F522;double F523;double F542;double F543;double XNO2;double AINV2;double TEMP1;double TEMP;double D2201;double D2211;double D3210;double D3222;double D4410;double D4422;double D5220;double D5232;double D5421;double D5433;double XLAMO;double BFACT;double G200;double G300;double F311;double F330;double DEL1;double DEL2;double DEL3;double FASX2;double FASX4;double FASX6;double XFACT;double XLI;double XNI;int IRET;int IRETN;double FT;double XNDOT;double XNDDT;double OMGDT;double XOMI;double X2OMI;double X2LI;double XLDOT;double XL;double T;double ZM;double ZF;double SINZF;double F2;double F3;double SES;double SIS;double SLS;double SGHS;double SHS;double SEL;double SIL;double SLL;double SGHL;double SHL;double PE;double PINC;double PL;double PGH;double PH;double SINIQ;double COSIQ;double SINOK;double COSOK;double ALFDP;double BETDP;double DALF;double DBET;double XLS;double DLS;double SINIS;double COSIS;void DPINIT(ElementSet es,double EQSQ,double SINIQ,double COSIQ,double RTEQSQ,double AO,double COSQ2,double SINOMO,double COSOMO,double BSQ,double XLLDOT,double OMGDT,double XNODOT,double XNODP){this.OMGDT = OMGDT;this.SINIQ = SINIQ;this.COSIQ = COSIQ;THGR=THETAG(es.DS50);EQ = es.EO;XNQ = XNODP;AQNV = 1.f/AO;XQNCL = es.XINCL;XMAO=es.XMO;XPIDOT=OMGDT+XNODOT;SINQ = Math.sin(es.XNODEO);COSQ = Math.cos(es.XNODEO);OMEGAQ = es.OMEGAO;// INITIALIZE LUNAR SOLAR TERMSDAY=es.DS50+18261.5;if (DAY != PREEP){PREEP = DAY;XNODCE=4.5236020-9.2422029E-4*DAY;STEM=Math.sin(XNODCE);CTEM=Math.cos(XNODCE);ZCOSIL=.91375164f-.03568096f*CTEM;ZSINIL=Math.sqrt(1.-ZCOSIL*ZCOSIL);ZSINHL= .089683511f*STEM/ZSINIL;ZCOSHL=Math.sqrt(1.-ZSINHL*ZSINHL);C=(4.7199672+.22997150*DAY);GAM=(5.8351514+.0019443680*DAY);ZMOL = FMOD2P(C-GAM);ZX= .39785416f*STEM/ZSINIL;ZY= ZCOSHL*CTEM+0.91744867f*ZSINHL*STEM;ZX=ACTAN(ZX,ZY);ZX=(GAM+ZX-XNODCE);ZCOSGL=Math.cos(ZX);ZSINGL=Math.sin(ZX);ZMOS=(6.2565837+.017201977*DAY);ZMOS=FMOD2P(ZMOS);}// DO SOLAR TERMSLS = 0;SAVTSN=1.e20;ZCOSG=ZCOSGS;ZSING=ZSINGS;ZCOSI=ZCOSIS;ZSINI=ZSINIS;ZCOSH=COSQ;ZSINH=SINQ;CC=C1SS;ZN=ZNS;ZE=ZES;ZMO=ZMOS;XNOI=1.f/XNQ;LS = 30;while (true){A1=ZCOSG*ZCOSH+ZSING*ZCOSI*ZSINH;A3=-ZSING*ZCOSH+ZCOSG*ZCOSI*ZSINH;A7=-ZCOSG*ZSINH+ZSING*ZCOSI*ZCOSH;A8=ZSING*ZSINI;A9=ZSING*ZSINH+ZCOSG*ZCOSI*ZCOSH;A10=ZCOSG*ZSINI;A2= COSIQ*A7+ SINIQ*A8;A4= COSIQ*A9+ SINIQ*A10;A5=- SINIQ*A7+ COSIQ*A8;A6=- SINIQ*A9+ COSIQ*A10;X1=A1*COSOMO+A2*SINOMO;X2=A3*COSOMO+A4*SINOMO;X3=-A1*SINOMO+A2*COSOMO;X4=-A3*SINOMO+A4*COSOMO;X5=A5*SINOMO;X6=A6*SINOMO;X7=A5*COSOMO;X8=A6*COSOMO;Z31=12.f*X1*X1-3.f*X3*X3;Z32=24.f*X1*X2-6.f*X3*X4;Z33=12.f*X2*X2-3.f*X4*X4;Z1=3.f*(A1*A1+A2*A2)+Z31*EQSQ;Z2=6.f*(A1*A3+A2*A4)+Z32*EQSQ;Z3=3.f*(A3*A3+A4*A4)+Z33*EQSQ;Z11=-6.f*A1*A5+EQSQ *(-24.f*X1*X7-6.f*X3*X5);Z12=-6.f*(A1*A6+A3*A5)+EQSQ *(-24.f*(X2*X7+X1*X8)-6.f*(X3*X6+X4*X5));Z13=-6.f*A3*A6+EQSQ *(-24.f*X2*X8-6.f*X4*X6);Z21=6.f*A2*A5+EQSQ *(24.f*X1*X5-6.f*X3*X7);Z22=6.f*(A4*A5+A2*A6)+EQSQ *(24.f*(X2*X5+X1*X6)-6.f*(X4*X7+X3*X8));Z23=6.f*A4*A6+EQSQ *(24.f*X2*X6-6.f*X4*X8);Z1=Z1+Z1+BSQ*Z31;Z2=Z2+Z2+BSQ*Z32;Z3=Z3+Z3+BSQ*Z33;S3=CC*XNOI;S2=-.5f*S3/RTEQSQ;S4=S3*RTEQSQ;S1=-15.f*EQ*S4;S5=X1*X3+X2*X4;S6=X2*X3+X1*X4;S7=X2*X4-X1*X3;SE=S1*ZN*S5;SI=S2*ZN*(Z11+Z13);SL=-ZN*S3*(Z1+Z3-14.f-6.f*EQSQ);SGH=S4*ZN*(Z31+Z33-6.f);SH=-ZN*S2*(Z21+Z23);if (XQNCL < 5.2359877E-2) SH=0.0f;EE2=2.f*S1*S6;E3=2.f*S1*S7;XI2=2.f*S2*Z12;XI3=2.f*S2*(Z13-Z11);XL2=-2.f*S3*Z2;XL3=-2.f*S3*(Z3-Z1);XL4=-2.f*S3*(-21.f-9.f*EQSQ)*ZE;XGH2=2.f*S4*Z32;XGH3=2.f*S4*(Z33-Z31);XGH4=-18.f*S4*ZE;XH2=-2.f*S2*Z22;XH3=-2.f*S2*(Z23-Z21);if (LS == 40) break;// DO LUNAR TERMSSSE = SE;SSI=SI;SSL=SL;SSH=SH/SINIQ;SSG=SGH-COSIQ*SSH;SE2=EE2;SI2=XI2;SL2=XL2;SGH2=XGH2;SH2=XH2;SE3=E3;SI3=XI3;SL3=XL3;SGH3=XGH3;SH3=XH3;SL4=XL4;SGH4=XGH4;LS=1;ZCOSG=ZCOSGL;ZSING=ZSINGL;ZCOSI=ZCOSIL;ZSINI=ZSINIL;ZCOSH=ZCOSHL*COSQ+ZSINHL*SINQ;ZSINH=SINQ*ZCOSHL-COSQ*ZSINHL;ZN=ZNL;CC=C1L;ZE=ZEL;ZMO=ZMOL;LS = 40;}SSE = SSE+SE;SSI=SSI+SI;SSL=SSL+SL;SSG=SSG+SGH-COSIQ/SINIQ*SH;SSH=SSH+SH/SINIQ;// GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITSIRESFL=0;ISYNFL=0;if (XNQ >= .0052359877 || XNQ <= .0034906585){if (XNQ < 8.26E-3 || XNQ > 9.24E-3) return;if (EQ < 0.5) return;IRESFL =1;EOC=EQ*EQSQ;G201=-.306f-(EQ-.64f)*.440f;if (EQ <= .65){G211=3.616f-13.247f*EQ+16.290f*EQSQ;G310=-19.302f+117.390f*EQ-228.419f*EQSQ+156.591f*EOC;G322=-18.9068f+109.7927f*EQ-214.6334f*EQSQ+146.5816f*EOC;G410=-41.122f+242.694f*EQ-471.094f*EQSQ+313.953f*EOC;G422=-146.407f+841.880f*EQ-1629.014f*EQSQ+1083.435f*EOC;G520=-532.114f+3017.977f*EQ-5740f*EQSQ+3708.276f*EOC;}else{G211=-72.099f+331.819f*EQ-508.738f*EQSQ+266.724f*EOC;G310=-346.844f+1582.851f*EQ-2415.925f*EQSQ+1246.113f*EOC;G322=-342.585f+1554.908f*EQ-2366.899f*EQSQ+1215.972f*EOC;G410=-1052.797f+4758.686f*EQ-7193.992f*EQSQ+3651.957f*EOC;G422=-3581.69f+16178.11f*EQ-24462.77f*EQSQ+12422.52f*EOC;if (EQ <= .715){G520=1464.74f-4664.75f*EQ+3763.64f*EQSQ;}else{G520=-5149.66f+29936.92f*EQ-54087.36f*EQSQ+31324.56f*EOC;}}if (EQ < .7){G533=-919.2277f+4988.61f*EQ-9064.77f*EQSQ+5542.21f*EOC;G521 = -822.71072f+4568.6173f*EQ-8491.4146f*EQSQ+5337.524f*EOC;G532 = -853.666f+4690.25f*EQ-8624.77f*EQSQ+5341.4f*EOC;}else{G533=-37995.78f+161616.52f*EQ-229838.2f*EQSQ+109377.94f*EOC;G521 = -51752.104f+218913.95f*EQ-309468.16f*EQSQ+146349.42f*EOC;G532 = -40023.88f+170470.89f*EQ-242699.48f*EQSQ+115605.82f*EOC;}SINI2=SINIQ*SINIQ;F220=.75f*(1.f+2.f*COSIQ+COSQ2);F221=1.5f*SINI2;F321=1.875f*SINIQ*(1.f-2.f*COSIQ-3.f*COSQ2);F322=-1.875f*SINIQ*(1.f+2.f*COSIQ-3.f*COSQ2);F441=35.f*SINI2*F220;F442=39.3750f*SINI2*SINI2;F522=9.84375f*SINIQ*(SINI2*(1.f-2.f*COSIQ-5.f*COSQ2)+.33333333f*(-2.f+4.f*COSIQ+6.f*COSQ2));F523 = SINIQ*(4.92187512f*SINI2*(-2.f-4.f*COSIQ+10.f*COSQ2)+6.56250012f*(1.f+2.f*COSIQ-3.f*COSQ2));F542 = 29.53125f*SINIQ*(2.f-8.f*COSIQ+COSQ2*(-12.f+8.f*COSIQ+10.f*COSQ2));F543=29.53125f*SINIQ*(-2.f-8.f*COSIQ+COSQ2*(12.f+8.f*COSIQ-10.f*COSQ2));XNO2=XNQ*XNQ;AINV2=AQNV*AQNV;TEMP1 = 3.f*XNO2*AINV2;TEMP = TEMP1*ROOT22;D2201 = TEMP*F220*G201;D2211 = TEMP*F221*G211;TEMP1 = TEMP1*AQNV;TEMP = TEMP1*ROOT32;D3210 = TEMP*F321*G310;D3222 = TEMP*F322*G322;TEMP1 = TEMP1*AQNV;TEMP = 2.f*TEMP1*ROOT44;D4410 = TEMP*F441*G410;D4422 = TEMP*F442*G422;TEMP1 = TEMP1*AQNV;TEMP = TEMP1*ROOT52;D5220 = TEMP*F522*G520;D5232 = TEMP*F523*G532;TEMP = 2.f*TEMP1*ROOT54;D5421 = TEMP*F542*G521;D5433 = TEMP*F543*G533;XLAMO = (XMAO+2.0f*es.XNODEO-2.0*THGR);BFACT = XLLDOT+XNODOT+XNODOT-THDT-THDT;BFACT=BFACT+SSL+SSH+SSH;}else{// SYNCHRONOUS RESONANCE TERMS INITIALIZATIONIRESFL=1;ISYNFL=1;G200=1.0f+EQSQ*(-2.5f+.8125f*EQSQ);G310=1.0f+2.0f*EQSQ;G300=1.0f+EQSQ*(-6.0f+6.60937f*EQSQ);F220=.75f*(1.f+COSIQ)*(1.f+COSIQ);F311=.9375f*SINIQ*SINIQ*(1.f+3.f*COSIQ)-.75f*(1.f+COSIQ);F330=1.f+COSIQ;F330=1.875f*F330*F330*F330;DEL1=3.f*XNQ*XNQ*AQNV*AQNV;DEL2=2.f*DEL1*F220*G200*Q22;DEL3=3.f*DEL1*F330*G300*Q33*AQNV;DEL1=DEL1*F311*G310*Q31*AQNV;FASX2=.13130908f;FASX4=2.8843198f;FASX6=.37448087f;XLAMO=(XMAO+es.XNODEO+es.OMEGAO-THGR);BFACT = XLLDOT+XPIDOT-THDT;BFACT=BFACT+SSL+SSG+SSH;}XFACT=BFACT-XNQ;// INITIALIZE INTEGRATORXLI=XLAMO;XNI=XNQ;ATIME=0.0;STEPP=720.0;STEPN=-720.0;STEP2 = 259200.0;}void DPSEC(ElementSet es,DoubleRef XLL,DoubleRef OMGASM,DoubleRef XNODES,DoubleRef EM,DoubleRef XINC,DoubleRef XN,double T){this.T = T;XLL.value=XLL.value+SSL*T;OMGASM.value=OMGASM.value+SSG*T;XNODES.value=XNODES.value+SSH*T;EM.value=es.EO+SSE*T;XINC.value=es.XINCL+SSI*T;if (XINC.value < 0.){XINC.value = -XINC.value;XNODES.value = XNODES.value + PI;OMGASM.value = OMGASM.value - PI;}if (IRESFL == 0) return;while (true){boolean checkIRETN = true;if ((ATIME == 0.0) ||(T >= 0.0 && ATIME < 0.0) ||(T < 0.0 && ATIME >= 0.0)){// EPOCH RESTARTif (T < 0.0){DELT = STEPN;}else{DELT = STEPP;}ATIME = 0.0;XNI=XNQ;XLI=XLAMO;}else{if (Math.abs(T) < Math.abs(ATIME)){DELT=STEPP;if (T >= 0.0) DELT = STEPN;IRET = 100;IRETN = 165;checkIRETN = false;}else{DELT=STEPN;if (T > 0.0) DELT = STEPP;}}while (true){if (checkIRETN){if (Math.abs(T-ATIME) >= STEPP){IRET = 125;IRETN = 165;}else{FT = (T-ATIME);IRETN = 140;}}checkIRETN = true;// DOT TERMS CALCULATEDif (ISYNFL != 0){XNDOT=DEL1*Math.sin(XLI-FASX2)+DEL2*Math.sin(2.f*(XLI-FASX4))+DEL3*Math.sin(3.f*(XLI-FASX6));XNDDT = DEL1*Math.cos(XLI-FASX2)+2.f*DEL2*Math.cos(2.f*(XLI-FASX4))+3.f*DEL3*Math.cos(3.f*(XLI-FASX6));}else{XOMI = (OMEGAQ+OMGDT*ATIME);X2OMI = XOMI+XOMI;X2LI = XLI+XLI;XNDOT = D2201*Math.sin(X2OMI+XLI-G22)+D2211*Math.sin(XLI-G22)+D3210*Math.sin(XOMI+XLI-G32)+D3222*Math.sin(-XOMI+XLI-G32)+D4410*Math.sin(X2OMI+X2LI-G44)+D4422*Math.sin(X2LI-G44)+D5220*Math.sin(XOMI+XLI-G52)+D5232*Math.sin(-XOMI+XLI-G52)+D5421*Math.sin(XOMI+X2LI-G54)+D5433*Math.sin(-XOMI+X2LI-G54);XNDDT = (D2201*Math.cos(X2OMI+XLI-G22)+D2211*Math.cos(XLI-G22)+D3210*Math.cos(XOMI+XLI-G32)+D3222*Math.cos(-XOMI+XLI-G32)+D5220*Math.cos(XOMI+XLI-G52)+D5232*Math.cos(-XOMI+XLI-G52)+2.*(D4410*Math.cos(X2OMI+X2LI-G44)+D4422*Math.cos(X2LI-G44)+D5421*Math.cos(XOMI+X2LI-G54)+D5433*Math.cos(-XOMI+X2LI-G54)));}XLDOT=XNI+XFACT;XNDDT = XNDDT*XLDOT;if (IRETN == 140){XN.value = XNI+XNDOT*FT+XNDDT*FT*FT*0.5f;XL = XLI+XLDOT*FT+XNDOT*FT*FT*0.5f;TEMP = (-XNODES.value+THGR+T*THDT);XLL.value = XL-OMGASM.value+TEMP;if (ISYNFL == 0) XLL.value = XL+TEMP+TEMP;return;}// INTEGRATORXLI = (XLI+XLDOT*DELT+XNDOT*STEP2);XNI = (XNI+XNDOT*DELT+XNDDT*STEP2);ATIME=ATIME+DELT;if (IRET == 100){break;}}}}// ENTRANCES FOR LUNAR-SOLAR PERIODICSvoid DPPER(DoubleRef EM,DoubleRef XINC,DoubleRef OMGASM,DoubleRef XNODES,DoubleRef XLL){SINIS = Math.sin(XINC.value);COSIS = Math.cos(XINC.value);if (Math.abs(SAVTSN-T) >= 30.0){SAVTSN=T;ZM=ZMOS+ZNS*T;ZF=ZM+2.f*ZES*Math.sin(ZM);SINZF=Math.sin(ZF);F2=.5f*SINZF*SINZF-.25f;F3=-.5f*SINZF*Math.cos(ZF);SES=SE2*F2+SE3*F3;SIS=SI2*F2+SI3*F3;SLS=SL2*F2+SL3*F3+SL4*SINZF;SGHS=SGH2*F2+SGH3*F3+SGH4*SINZF;SHS=SH2*F2+SH3*F3;ZM=ZMOL+ZNL*T;ZF=ZM+2.f*ZEL*Math.sin(ZM);SINZF=Math.sin(ZF);F2=.5f*SINZF*SINZF-.25f;F3=-.5f*SINZF*Math.cos(ZF);SEL=EE2*F2+E3*F3;SIL=XI2*F2+XI3*F3;SLL=XL2*F2+XL3*F3+XL4*SINZF;SGHL=XGH2*F2+XGH3*F3+XGH4*SINZF;SHL=XH2*F2+XH3*F3;PE=SES+SEL;PINC=SIS+SIL;PL=SLS+SLL;}PGH=SGHS+SGHL;PH=SHS+SHL;XINC.value = XINC.value+PINC;EM.value = EM.value+PE;if (XQNCL >= .2){// APPLY PERIODICS DIRECTLYPH=PH/SINIQ;PGH=PGH-COSIQ*PH;OMGASM.value=OMGASM.value+PGH;XNODES.value=XNODES.value+PH;XLL.value = XLL.value+PL;}else{// APPLY PERIODICS WITH LYDDANE MODIFICATIONSINOK=Math.sin(XNODES.value);COSOK=Math.cos(XNODES.value);ALFDP=SINIS*SINOK;BETDP=SINIS*COSOK;DALF=PH*COSOK+PINC*COSIS*SINOK;DBET=-PH*SINOK+PINC*COSIS*COSOK;ALFDP=ALFDP+DALF;BETDP=BETDP+DBET;XLS = XLL.value+OMGASM.value+COSIS*XNODES.value;DLS=PL+PGH-PINC*XNODES.value*SINIS;XLS=XLS+DLS;XNODES.value=ACTAN(ALFDP,BETDP);XLL.value = XLL.value+PL;OMGASM.value = XLS-XLL.value-Math.cos(XINC.value)*XNODES.value;}}}
class Earth{int latRes = 20;int lonRes = 20;float rad = (float)XKMPER;Tuple3 [][] pts;public Earth(){pts = new Tuple3[lonRes][latRes+1];float lat = 0;float lon = 0;for (int i = 0; i < lonRes; i++){lon = (float)PI + (float)(i*2*PI)/lonRes;for (int j = 0; j < latRes+1; j++){lat = (float)(j*(PI/latRes));pts[i][j] = new Tuple3();pts[i][j].x = rad*cos(lon)*sin(lat);pts[i][j].y = rad*sin(lon)*sin(lat);pts[i][j].z = rad*cos(lat);pts[i][j].uv = (float)i/lonRes;pts[i][j].uw = (float)j/latRes;}}}void draw(){textureMode(NORMALIZED);Tuple3 pt, pt2,pt3,pt4;pushMatrix();rotateZ((float)lst);((PGraphics3D)g).triangle.setCulling(true);//fast render loop for unshaded earth//fill(255,255,255,230);//if (false){tint(255,255,255);beginShape(TRIANGLE_STRIP);texture(img);for (int i = 0; i < latRes; i++){for (int j = 0; j < lonRes; j++){pt = pts[j][i+1];pt2 = pts[j][i];//tint(random(100,200),100,100);vertex(pt.x,pt.y,pt.z,pt.uv,pt.uw);vertex(pt2.x,pt2.y,pt2.z,pt2.uv,pt2.uw);}//close the looppt = pts[0][i];pt2 = pts[0][i+1];vertex(pt2.x,pt2.y,pt2.z,1,pt2.uw);vertex(pt.x,pt.y,pt.z,1,pt.uw);}endShape();/*}else{//slow render loop for night skyfor (int i = 0; i < latRes; i++){for (int j = 0; j < lonRes; j++){pt3 = pts[j][i];pt4 = pts[j][i+1];pt = pts[j+1][i];pt2 = pts[j+1][i+1];beginShape(TRIANGLES);if (pt.x > 0) texture(img);else texture(nimg);normal(pt.x,pt.y,pt.z);vertex(pt.x,pt.y,pt.z,pt.uv,pt.uw);normal(pt2.x,pt2.y,pt2.z);vertex(pt2.x,pt2.y,pt2.z,pt2.uv,pt2.uw);normal(pt3.x,pt3.y,pt3.z);vertex(pt3.x,pt3.y,pt3.z,pt3.uv,pt3.uw);normal(pt2.x,pt2.y,pt2.z);vertex(pt2.x,pt2.y,pt2.z,pt2.uv,pt2.uw);normal(pt3.x,pt3.y,pt3.z);vertex(pt3.x,pt3.y,pt3.z,pt3.uv,pt3.uw);normal(pt4.x,pt4.y,pt4.z);vertex(pt4.x,pt4.y,pt4.z,pt4.uv,pt4.uw);endShape();}}}*/popMatrix();((PGraphics3D)g).triangle.setCulling(false);if (showEquator == 1){fill(255,255,255,50);pushMatrix();ellipse(0,0,2*(35786 + erad),2*(35786 + erad));popMatrix();}}}class Tuple3{float x,y,z;float uv,uw;}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*//****/public class ElementSet {double XMO = 0.0;double XNODEO = 0.0;double OMEGAO = 0.0;double EO = 0.0;double XINCL = 0.0;double XNO = 0.0;double XNDT2O = 0.0;double XNDD6O = 0.0;double BSTAR = 0.0;double EPOCH = 0.0;double DS50 = 0.0;double julian = 0.0;int satNum;private final String line1;private final String line2;/*** Constructs an instance of this class from a NORAD two-line element set.*/public ElementSet(String[] tle) {this(tle[0],tle[1]);}/*** Constructs an instance of this class from a NORAD two-line element set.*/public ElementSet(String tleLine1,String tleLine2) {this.line1 = tleLine1;this.line2 = tleLine2;int IEXP =0;int IBEXP = 0;satNum = Integer.parseInt(line1.substring(2,7));//line1EPOCH = Double.parseDouble(line1.substring(18,32));int year = Integer.parseInt(line1.substring(18,20));double epoch = Double.parseDouble(line1.substring(20,32));julian = toJulian(year,epoch);XNDT2O = Double.parseDouble(line1.substring(34,43));XNDD6O = Double.parseDouble("." + line1.substring(45,50));IEXP = Integer.parseInt(line1.substring(50,52));BSTAR = Double.parseDouble("." + line1.substring(54,59));if (line1.substring(59,60).equals("+")){IBEXP = Integer.parseInt(line1.substring(60,61));}else{IBEXP = Integer.parseInt(line1.substring(59,61));}XINCL = Double.parseDouble(line2.substring(8,16));XNODEO = Double.parseDouble(line2.substring(17,25));EO = Double.parseDouble("." + line2.substring(26,33));OMEGAO = Double.parseDouble(line2.substring(34,42));XMO = Double.parseDouble(line2.substring(43,51));XNO = Double.parseDouble(line2.substring(52,63));XNDD6O=XNDD6O*Math.pow(10.,IEXP);XNODEO=XNODEO*DE2RA;OMEGAO=OMEGAO*DE2RA;XMO=XMO*DE2RA;XINCL=XINCL*DE2RA;double TEMP=TWOPI/XMNPDA/XMNPDA;XNO=XNO*TEMP*XMNPDA;XNDT2O=XNDT2O*TEMP;XNDD6O=XNDD6O*TEMP/XMNPDA;BSTAR=BSTAR*Math.pow(10.,IBEXP)/AE;DS50=computeDS50(EPOCH);}/*** @return Line 2 of the two-line element set.*/public String getLine1() {return line1;}/*** @return Line 2 of the two-line element set.*/public String getLine2() {return line2;}/*** @return true if the element set is for a deep-space orbit, false if it's* for a near-earth orbit.*/public boolean isDeep() {double A1=Math.pow((XKE/XNO),TOTHRD);double TEMP=1.5*CK2*(3.*Math.pow(Math.cos(XINCL),2)-1.)/Math.pow((1.-EO*EO),1.5);double DEL1=TEMP/(A1*A1);double AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)));double DELO=TEMP/(AO*AO);double XNODP=XNO/(1.+DELO);return ((TWOPI/XNODP/XMNPDA) >= .15625);}}
import processing.net.*;import javax.swing.*;Camera cam;PFont font;int numStars = 1500;Tuple [] stars = new Tuple[numStars];String lines[];ElementSet es;SGP4 propagator;Vector history = new Vector();//the amount the earth rotates in radians per time stepdouble lstStep = 0;int showEquator = 0;int Y,M,D;float mod = 1.25;Vector sats = new Vector();Earth earth = new Earth();int displayHelp = 0;double trailLength = .025;//stupid eyecandy and its first order lagfloat controlBlur = 0;float dControlBlur = 0;String [] months = {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec"};int pickMode = 0;int camMode = 0;PImage img,nimg;double timeSpan = 24;double timeStep = 3;Propagator current = null;int showLabels = 1;int geoTrack = 0;int realTime = 0;int showInfo = 1;int running = 1;float ts = 1; //in minutesdouble lst = 0;double JD = 0;int year = 0;double day = 0;int hour = 0;int minute = 0;int second = 0;float longitude =0;float time_zone = 0;double runStep = .00005;double startDate = 0;//this is astro, we work in GMTCalendar calendar = new GregorianCalendar();String onLoad = null;void setup(){String TITLE="0";int w,h;if (online){w = int(param("width"));if (w <= 0) {w = 550;}h = int(param("height"));if (h <= 0) {h = 550;}try{URL base = getDocumentBase();StringTokenizer st = new StringTokenizer(base.toString(),"?=");String a = st.nextToken();String par = st.nextToken();onLoad = st.nextToken();if (par.equals("TITLE")){onLoad = st.nextToken();}}catch(Exception e){println("No Input values");}}else{String lines[] = loadStrings("config.txt");if (lines != null){w = new Integer(lines[0]).intValue();h = new Integer(lines[1]).intValue();}else{w = 600;h = 600;}}size(w,h,P3D);JD = getCurrentJulianDate();startDate = JD;img = loadImage("earth.jpg");cam = new Camera(0,20000,0, 0,0,0);font = loadFont("writing.vlw");textFont(font);textMode(SCREEN);for (int i = 0; i < numStars; i++){stars[i] = new Tuple(random(-5000000,5000000),random(-5000000,5000000),random(-5000000,5000000));}if (onLoad != null ){if (onLoad.equals("Iridium")){loadFile("iridium.txt");showLabels = 0;trailLength = .0075;}else if (onLoad.equals("GPS")){loadFile("gps.txt");cam.cy=60000;cam.calcShotLength();}else if (onLoad.equals("Molniya")){loadFile("molniya.txt");cam.cy=80000;cam.calcShotLength();runStep = .0002;}else if (onLoad.equals("TDRS")){loadFile("tdrss.txt");cam.cy=90000;cam.calcShotLength();}}else{loadFile("recent.txt");}lstStep = 2*PI*(timeStep/1436.06818);}void draw(){background(0);if (camMode == 1 ){if (current != null){cam.tx = -current.currX;cam.ty = -current.currY;cam.tz = -current.currZ;cam.calcShotLength();}else{camMode =0;}}cam.feed();if (running == 1) JD = JD + runStep;if (geoTrack == 1 && running == 1 && realTime == 0) cam.rotZ(-(float)(2*PI*runStep));if (realTime == 1) JD = getCurrentJulianDate();lst = getSiderealTime(JD);stroke(150);beginShape(POINTS);for (int i = 0; i < numStars; i++){vertex(stars[i].x,stars[i].y,stars[i].z);}endShape();noStroke();fill(100,100,200);earth.draw();stroke(255);noFill();Propagator prop;pushMatrix();for (int i =0; i < sats.size(); i++){prop = (Propagator)sats.elementAt(i);prop.render();}popMatrix();//let'g get back to some screen coord goodnesscamera();stroke(128);fill(200,200,255,50);if (current != null){stroke(255,255,0);noFill();ellipse((int)current.msx,(int)current.msy,30,30);}//where are we in our time cache?//in julian daysif (mouseY > (height - 50) ) dControlBlur = 150;else dControlBlur = 50;//negative feedbackcontrolBlur += (dControlBlur - controlBlur)/10;double gaugeWidth = timeSpan/24;int rectWidth = width-100;double offset = JD-startDate;double percent = offset/gaugeWidth;stroke(100);if (percent > 1) {percent = 1;running = 0;}if (percent < 0){percent = 0;running = 0;}fill(100,100,120,controlBlur);rect(50,height-30,rectWidth,20);fill(200,200,200,controlBlur);int base = (int)(percent*(rectWidth-20));rect(50 + base,height-30,20,20);if (showInfo ==1){//legibility boxfill(0,0,200,50);rect(5,5,230,55);fill(255);if (realTime == 1) text("Realtime Julian Date: " + JD,10,15);else text("Displayed Julian Date: " + JD,10,15);//text(formatter.format( refCalendar.getTime()) + " GMT " + hour + ":" + minute + ":" + second ,10,25);text(D + " " + months[M-1] + " " + Y + " GMT " + hour + ":" + minute + ":" + second ,10,25);text(timeSpan + " hours of ephemeris",10,35);if (geoTrack == 1) text("Camera rotates with Earth",10,45);else text("Camera is inertial",10,45);text("Time Step: " + (float)(runStep*86400) + " seconds per tick",10,55);//have to do one set of render/calculates}if (pickMode == 1) pickMode = 2;if (displayHelp == 1){fill(0,0,0,200);noStroke();rect(0,0,width,height);fill(255,255,255);int ht = 50;int wd = width/2 - 100;int sprd = 15;text("1 - Load ASAT debris",wd,ht);ht += sprd;text("2 - Load ISS",wd,ht);ht += sprd;text("3 - Load Satellite Radio",wd,ht);ht += sprd;text("4 - Load GPS",wd,ht);ht += sprd;text("5 - Load Recent launches",wd,ht);ht += sprd;text("6 - Load Molniya orbits",wd,ht);ht += sprd;text("7 - Load Iridium",wd,ht);ht += sprd;text("8 - Load TDRS",wd,ht);ht += sprd;text("9 - Load Cubesats",wd,ht);ht += sprd;text("a - Load TLE from disk (Standalone version only)",wd,ht);ht += sprd;ht+= sprd;text("left mouse drag - Orbit",wd,ht);ht += sprd;text("right mouse drag - Zoom",wd,ht);ht += sprd;text("z + mouse drag - Zoom",wd,ht);ht += sprd;text("c - Center camera on Earth",wd,ht);ht += sprd;text("q - Remove all satellites",wd,ht);ht += sprd;text("p - Enter pick mode (click to select)",wd,ht);ht += sprd;text("u - Deselect all",wd,ht);ht += sprd;ht += sprd;text("g - Toggle inertial/rotate with earth",wd,ht);ht += sprd;text("up/down - Increase/decrease timestep",wd,ht);ht += sprd;text("left/right - move forward/backward in time",wd,ht);ht += sprd;text("r - Lock sim time to actual time",wd,ht);ht += sprd;text("space - Toggle running",wd,ht);ht += sprd;text("n - Set clock to current time",wd,ht);ht += sprd;ht+= sprd;text("l - Toggle labels",wd,ht);ht += sprd;text("d - Toggle infobox",wd,ht);ht += sprd;text("e - Toggle equator",wd,ht);ht += sprd;text("w/s - Increment/decrement trail length",wd,ht);ht += sprd;}}void mousePressed(){if (pickMode == 2){pick();pickMode = 0;}}void mouseDragged() {//quick hack so mac users with no second mouse button can zoomif (mouseButton == LEFT && key != 'z'){if (mouseY < height - 30){cam.upZ(-radians(mouseY - pmouseY)/2);cam.rotZ(radians(mouseX - pmouseX)/2);}else{realTime = 0;//how much time does pixel equaldouble gaugeWidth = timeSpan/24;int rectWidth = width-100;double pixelEquivalent = gaugeWidth/rectWidth;JD += pixelEquivalent*(mouseX-pmouseX);}}else if (mouseButton == RIGHT || key == 'z'){cam.modShot(1 + (mouseY - pmouseY)/50.f);}}void keyReleased(){key = ' ';}void keyPressed(){if (key != 'h') displayHelp = 0;if(key == ' ' ){realTime = 0;running = 1-running;println(running);}else if(key == 'n' ){JD = getCurrentJulianDate();}else if(key == '1' ){loadFile("ASAT.txt");showLabels = 0;}else if(key == '2' ){loadFile("ISS.txt");}else if(key == '3' ){loadFile("radio.txt");}else if(key == '4' ){loadFile("gps.txt");}else if(key == '5' ){loadFile("recent.txt");}else if(key == '6' ){loadFile("molniya.txt");}else if(key == '7' ){loadFile("iridium.txt");}else if(key == '8' ){loadFile("tdrss.txt");}else if(key == '9' ){loadFile("cubesat.txt");}else if(key == '0' ){loadFile("glonass.txt");}else if(key == 'd' ){showInfo = 1-showInfo;}else if(key == 'l' ){showLabels = 1-showLabels;}else if(key == 'r' ){realTime = 1-realTime;running = 0;}else if (key == 'p'){pickMode = 1;}else if (key == 'g'){geoTrack = 1- geoTrack;}else if (key == 'h'){displayHelp = 1- displayHelp;}else if (key == 'u'){if (current != null)current.selected = 0;current = null;camMode = 0;cam.tx = 0;cam.ty = 0;cam.tz = 0;cam.calcShotLength();}else if (keyCode == RIGHT){JD += runStep;}else if (keyCode == LEFT){JD -= runStep;}else if (keyCode == UP){runStep *= 1.1;}else if (keyCode == DOWN){runStep /= 1.1;}else if (key == 'c'){if (camMode ==1){camMode = 0;cam.tx = 0;cam.ty = 0;cam.tz = 0;cam.calcShotLength();}else if (current != null){camMode =1;}}else if (key == 'a'){appendFile();}else if (key == 'q'){sats = new Vector();current = null;}else if (key == 'w'){trailLength *= 1.1;}else if (key == 's'){trailLength /= 1.1;}else if (key == 'e'){showEquator = 1- showEquator;}try{if (key == 'S'){img = loadImage("earth_summer.jpg");}else if (key == 'W'){img = loadImage("earth_winter.jpg");}else if (key == 'U'){img = loadImage("earth_ultra.jpg");}}//image may not existcatch(Exception e){println("Could not load image");};}void pick(){Propagator prop;float minDist = 200000;float distance;for (int i =0; i < sats.size(); i++){prop = (Propagator)sats.elementAt(i);distance = prop.pick(mouseX,mouseY);if (distance < minDist){minDist = distance;if (current != null){current.selected = 0;}prop.selected = 1 ;current = prop;}}//camMode = 1;}void appendFile(){SwingUtilities.invokeLater(new Runnable() {public void run() {try {JFileChooser fc = new JFileChooser();int returnVal = fc.showOpenDialog(null);if (returnVal == JFileChooser.APPROVE_OPTION) {File selFile = fc.getSelectedFile();if (selFile != null){String outName = selFile.getAbsolutePath();loadFile(outName);}}} catch (Exception e) {e.printStackTrace();}}});}void loadFile(String fileName){//open the file with all the gps data in itlines = loadStrings(fileName);parseStrings();}void parseStrings(){try{Propagator prop;ElementSet es;//load all the satellites in the filefor (int i=0; i < lines.length; i+=3) {//watch out for blank lineswhile(lines[i] == "") i++;//use the celestrak (NORAD TLE) file formatprintln("Generating:" + lines[i]);es = new ElementSet(lines[i+1], lines[i+2]);if (es.isDeep()){prop = new SDP4(lines[i],new ElementSet(lines[i+1], lines[i+2]));}else{prop = new SGP4(lines[i],new ElementSet(lines[i+1], lines[i+2]));}//propagate the satellite and store within the propagatorprop.generate(startDate,timeSpan,timeStep);//perform a sanity check for de-orbited satellites.OrbitalState a,b;a = (OrbitalState)prop.history.elementAt(10);b = (OrbitalState)prop.history.elementAt(11);double dx = a.X-b.X;double dy = a.Y-b.Y;double dz = a.Z-b.Z;//sats don't move this fastif (Math.sqrt(dx*dx + dy*dy + dz*dz) > 10000){println("Excluding " + lines[i] + ": orbital motion incorrect." );}else{sats.add(prop);}}}catch(Exception e){println(e);}}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*/class DoubleRef {double value = 0;/*** Constructs an instance of this class which has the given value.*/DoubleRef(double value) {this.value = value;}}public static double ACTAN(double SINX, double COSX) {double ACTAN = 0.0f;if (COSX != 0. ){if (COSX <= 0. ){ACTAN=PI;}else{if (SINX == 0. ){return ACTAN;}if (SINX <= 0. ){ACTAN=TWOPI;}}}else{if (SINX == 0. ){return ACTAN;}if (SINX <= 0. ){ACTAN=X3PIO2;return ACTAN;}else{ACTAN=PIO2;return ACTAN;}}double TEMP=SINX/COSX;ACTAN=ACTAN+Math.atan(TEMP);return ACTAN;}public static double FMOD2P(double X) {double FMOD2P = 0.0;FMOD2P=X;int I=(int)(FMOD2P/TWOPI);FMOD2P=FMOD2P-I*TWOPI;if (FMOD2P < 0) FMOD2P=FMOD2P+TWOPI;return FMOD2P;}public static double computeDS50(double EPOCH) {double YR=(EPOCH+2.e-7)*1.e-3;int JY=(int)YR;YR=JY;double D=EPOCH-YR*1.e3;if (JY < 10) JY=JY+80;int N=(JY-69)/4;if (JY < 70) N=(JY-72)/4;return 7305. + 365.*(JY-70) + N + D;}public static double THETAG(double DS50) {double THETA=1.72944494 + 6.3003880987*DS50;double THETAG = THETA - ((int)(THETA/TWOPI))*TWOPI;if (THETAG < 0.0){THETAG = THETAG + TWOPI;}return THETAG;}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*/import java.text.DecimalFormat;/*** Represents orbital state of the spacecraft at a point in time.*/public final class OrbitalState {/*** Timestamp of the orbital state, in seconds since the epoch of the known* orbital state.*/public double julianTime = 0.0;/*** Position of the spacecraft at the point in time.*/public double X = 0.0;public double Y = 0.0;public double Z = 0.0;/*** Velocity of the spacecraft at the point in time.*/public double XDOT = 0.0;public double YDOT = 0.0;public double ZDOT = 0.0;/*** Constructs an instance of this class.*/public OrbitalState() {}public String toString(){return X + " : " + Y + ":" + Z;}}
abstract class Propagator{String title;Vector history = new Vector();double startTime;double stepSize;protected ElementSet es;float currX=0,currY=0,currZ=0;int selected = 0;float msx=0,msy=0;float temp;float az=0,lastAz=0;//used for mouse pickingfloat pick(int sx, int sy){return sqrt((msx-sx)*(msx-sx) + (msy-sy)*(msy-sy));}void render(){noFill();beginShape();stroke(255,255,255,100);OrbitalState os = null, nextOs = new OrbitalState(), currOs = new OrbitalState();double off;double perc = 0.0;//head interpolation point for trail and path generationint start = 0;boolean brk = false;for (int i =0; i < history.size()-1 && brk == false; i++){os = (OrbitalState)history.elementAt(i);off = JD - os.julianTime;if (off > 0 && off < trailLength){stroke(255,255,255,150 - (int)(150*off/trailLength));vertex((float)os.X,(float)os.Y,(float)os.Z);}if (off > 0 && off < stepSize){currOs = os;//these are used to interpolateperc = off/stepSize;nextOs=(OrbitalState)history.elementAt(i+1);currX = (float)(currOs.X + perc*(nextOs.X - currOs.X));currY = (float)(currOs.Y + perc*(nextOs.Y - currOs.Y));currZ = (float)(currOs.Z + perc*(nextOs.Z - currOs.Z));vertex(currX,currY,currZ);brk = true;start = i;}}endShape();//render groundPathif (selected == 1){double gx,gy,gz;double ngx;beginShape();stroke(255,255,0);double ht;double inv;//step backwardsht = sqrt(currX*currX + currY*currY + currZ*currZ);vertex(currX,currY,currZ);vertex((float)(erad*currX/ht),(float)(erad*currY/ht),(float)(erad*currZ/ht));for (int i = start; i >0; i--){os = (OrbitalState)history.elementAt(i);inv = (start - i + perc);//get height to normalize pointing vectorht = Math.sqrt(os.X*os.X + os.Y*os.Y + os.Z*os.Z);gx = os.X;gy = os.Y;gz = os.Z;//unwrap the groundpath based on the earth's past rotation (averaged)ngx = gx*Math.cos(inv*lstStep) - gy*Math.sin(inv*lstStep);gy = gx*Math.sin(inv*lstStep) + gy*Math.cos(inv*lstStep);ht = Math.sqrt(ngx*ngx + gy*gy + gz*gz);ngx /= ht;gy /= ht;gz /= ht;ngx*= 1.005*erad;gy *= 1.005*erad;gz *= 1.005*erad;vertex((float)ngx,(float)gy,(float)gz);}endShape();/*//draw the groundpath windowfill(255,255,255,255);int wi = width/2;int hi = (img.height*wi)/img.width;pushMatrix();camera();stroke(255);translate(width - wi - 10, 10);beginShape(QUADS);texture(img);vertex(0,0,0,0);vertex(wi,0,1,0);vertex(wi,hi,1,1);vertex(0,hi,0,1);endShape();popMatrix();*/}//if we found the satelliteif (currX != 0){pushMatrix();translate(currX,currY,currZ);fill(200,200,255,100);stroke(200);box(200);fill(255);if (selected == 1 || pickMode == 1){msx = screenX(0,0,0);msy = screenY(0,0,0);lastAz = az;az = atan2(currX,currY);}if (showLabels == 1 || selected == 1){//text(title + ":" + es.satNum,screenX(0,0,0),screenY(0,0,0));text(title,screenX(0,0,0),screenY(0,0,0));}popMatrix();}if (false){beginShape();stroke(255);vertex(currX/erad,currY/erad,currZ/erad);vertex(currX,currY,currZ);endShape();}}OrbitalState propagate(double TSINCE){return new OrbitalState();}//generate this satellites ephemerides for the appointed time spanvoid generate(double fromJD, double span, double resolution){double offset = fromJD - startTime + mod;//convert offset from julian days to minutesoffset *= 24*60;//convert span from hours to minutesspan*=60;stepSize = resolution/(24*60);for (int j =0; j < span/resolution; j++){history.add(propagate(offset + j*resolution));}}}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*/public class SDP4 extends Propagator{private double A1;private double COSIO;private double THETA2;private double X3THM1;private double EOSQ;private double BETAO2;private double BETAO;private double DEL1;private double DELO;private double AO;private double XNODP;private double AODP;private double S4;private double QOMS24;private double PERIGE;private double PINVSQ;private double SING;private double COSG;private double TSI;private double ETA;private double ETASQ;private double EETA;private double PSISQ;private double COEF;private double COEF1;private double c1;private double c2;private double C4;private double SINIO;private double A3OVK2;private double X1MTH2;private double THETA4;private double TEMP1;private double TEMP2;private double TEMP3;private double XMDOT;private double X1M5TH;private double OMGDOT;private double XHDOT1;private double XNODOT;private double XNODCF;private double T2COF;private double XLCOF;private double AYCOF;private double X7THM1;private DoubleRef XMDF = new DoubleRef(0.0);private DoubleRef OMGADF = new DoubleRef(0.0);private double XNODDF;private double TSQ;private DoubleRef XNODE = new DoubleRef(0.0);private double TEMPA;private double TEMPE;private double TEMPL;private DoubleRef XN = new DoubleRef(0.0);private DoubleRef EM = new DoubleRef(0.0);private DoubleRef XINC = new DoubleRef(0.0);private double A;private DoubleRef E = new DoubleRef(0.0);private DoubleRef XMAM = new DoubleRef(0.0);private double XL;private double BETA;private double AXN;private double TEMP;private double XLL;private double AYNL;private double XLT;private double AYN;private double CAPU;private double SINEPW;private double COSEPW;private double TEMP4;private double TEMP5;private double TEMP6;private double EPW;private double ECOSE;private double ESINE;private double ELSQ;private double PL;private double R;private double RDOT;private double RFDOT;private double BETAL;private double COSU;private double SINU;private double U;private double COS2U;private double SIN2U;private double RK;private double UK;private double XNODEK;private double XINCK;private double RDOTK;private double RFDOTK;private double SINUK;private double COSUK;private double SINIK;private double COSIK;private double SINNOK;private double COSNOK;private double XMX;private double XMY;private double UX;private double UY;private double UZ;private double VX;private double VY;private double VZ;private DEEP deep;/*** Constructs an instance of this class.*/public SDP4(String title,ElementSet es) {this.title = title;this.es = es;startTime = es.julian;initialize();}private void initialize() {// RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)// FROM INPUT ELEMENTSA1=Math.pow((XKE/es.XNO),TOTHRD);COSIO=Math.cos(es.XINCL);THETA2=COSIO*COSIO;X3THM1=3.*THETA2-1.;EOSQ=es.EO*es.EO;BETAO2=1.-EOSQ;BETAO=Math.sqrt(BETAO2);DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2);AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)));DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2);XNODP=es.XNO/(1.+DELO);AODP=AO/(1.-DELO);// INITIALIZATION// FOR PERIGEE BELOW 156 KM, THE VALUES OF// S AND QOMS2T ARE ALTEREDS4=S;QOMS24=QOMS2T;PERIGE=(AODP*(1.-es.EO)-AE)*XKMPER;if (PERIGE < 156.){S4=PERIGE-78.;if (PERIGE <= 98.){S4=20.;}QOMS24=Math.pow(((120.-S4)*AE/XKMPER),4);S4=S4/XKMPER+AE;}PINVSQ=1./(AODP*AODP*BETAO2*BETAO2);SING=Math.sin(es.OMEGAO);COSG=Math.cos(es.OMEGAO);TSI=1./(AODP-S4);ETA=AODP*es.EO*TSI;ETASQ=ETA*ETA;EETA=es.EO*ETA;PSISQ=Math.abs(1.-ETASQ);COEF=QOMS24*Math.pow(TSI,4);COEF1=COEF/Math.pow(PSISQ,3.5);c2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)));c1=es.BSTAR*c2;SINIO=Math.sin(es.XINCL);A3OVK2=-XJ3/CK2*Math.pow(AE,3);X1MTH2=1.-THETA2;C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*(2.+.5*ETASQ)+es.EO*(.5+2.*ETASQ)-2.*CK2*TSI/(AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*(1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*(1.+ETASQ))*Math.cos(2.*es.OMEGAO)));THETA4=THETA2*THETA2;TEMP1=3.*CK2*PINVSQ*XNODP;TEMP2=TEMP1*CK2*PINVSQ;TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP;XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*(13.-78.*THETA2+137.*THETA4);X1M5TH=1.-5.*THETA2;OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4);XHDOT1=-TEMP1*COSIO;XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-7.*THETA2))*COSIO;XNODCF=3.5*BETAO2*XHDOT1*c1;T2COF=1.5*c1;XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO);AYCOF=.25*A3OVK2*SINIO;X7THM1=7.*THETA2-1.;deep = new DEEP();deep.DPINIT(es,EOSQ,SINIO,COSIO,BETAO,AODP,THETA2,SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP);}public OrbitalState propagate(double TSINCE) {// UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAGXMDF.value=es.XMO+XMDOT*TSINCE;OMGADF.value=es.OMEGAO+OMGDOT*TSINCE;XNODDF=es.XNODEO+XNODOT*TSINCE;TSQ=TSINCE*TSINCE;XNODE.value=XNODDF+XNODCF*TSQ;TEMPA=1.-c1*TSINCE;TEMPE=es.BSTAR*C4*TSINCE;TEMPL=T2COF*TSQ;XN.value=XNODP;deep.DPSEC(es,XMDF,OMGADF,XNODE,EM,XINC,XN,TSINCE);A=Math.pow((XKE/XN.value),TOTHRD)*Math.pow(TEMPA,2);E.value=EM.value-TEMPE;XMAM.value=XMDF.value+XNODP*TEMPL;deep.DPPER(E,XINC,OMGADF,XNODE,XMAM);XL=XMAM.value+OMGADF.value+XNODE.value;BETA=Math.sqrt(1.-E.value*E.value);XN.value=XKE/Math.pow(A,1.5);// LONG PERIOD PERIODICSAXN=E.value*Math.cos(OMGADF.value);TEMP=1./(A*BETA*BETA);XLL=TEMP*XLCOF*AXN;AYNL=TEMP*AYCOF;XLT=XL+XLL;AYN=E.value*Math.sin(OMGADF.value)+AYNL;// SOLVE KEPLERS EQUATIONCAPU=FMOD2P(XLT-XNODE.value);TEMP2=CAPU;for ( int I = 1 ; I <= 10 ; I++ ){SINEPW=Math.sin(TEMP2);COSEPW=Math.cos(TEMP2);TEMP3=AXN*SINEPW;TEMP4=AYN*COSEPW;TEMP5=AXN*COSEPW;TEMP6=AYN*SINEPW;EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2;if (Math.abs(EPW-TEMP2) <= E6A) break;TEMP2=EPW;}// SHORT PERIOD PRELIMINARY QUANTITIESECOSE=TEMP5+TEMP6;ESINE=TEMP3-TEMP4;ELSQ=AXN*AXN+AYN*AYN;TEMP=1.-ELSQ;PL=A*TEMP;R=A*(1.-ECOSE);TEMP1=1./R;RDOT=XKE*Math.sqrt(A)*ESINE*TEMP1;RFDOT=XKE*Math.sqrt(PL)*TEMP1;TEMP2=A*TEMP1;BETAL=Math.sqrt(TEMP);TEMP3=1./(1.+BETAL);COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3);SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3);U=ACTAN(SINU,COSU);SIN2U=2.*SINU*COSU;COS2U=2.*COSU*COSU-1.;TEMP=1./PL;TEMP1=CK2*TEMP;TEMP2=TEMP1*TEMP;// UPDATE FOR SHORT PERIODICSRK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U;UK=U-.25*TEMP2*X7THM1*SIN2U;XNODEK=XNODE.value+1.5*TEMP2*COSIO*SIN2U;XINCK=XINC.value+1.5*TEMP2*COSIO*SINIO*COS2U;RDOTK=RDOT-XN.value*TEMP1*X1MTH2*SIN2U;RFDOTK=RFDOT+XN.value*TEMP1*(X1MTH2*COS2U+1.5*X3THM1);// ORIENTATION VECTORSSINUK=Math.sin(UK);COSUK=Math.cos(UK);SINIK=Math.sin(XINCK);COSIK=Math.cos(XINCK);SINNOK=Math.sin(XNODEK);COSNOK=Math.cos(XNODEK);XMX=-SINNOK*COSIK;XMY=COSNOK*COSIK;UX=XMX*SINUK+COSNOK*COSUK;UY=XMY*SINUK+SINNOK*COSUK;UZ=SINIK*SINUK;VX=XMX*COSUK-COSNOK*SINUK;VY=XMY*COSUK-SINNOK*SINUK;VZ=SINIK*COSUK;OrbitalState sv = new OrbitalState();//convert minutes elapsed into fractional days.double deltaTime = TSINCE*60/86400 - mod ;sv.julianTime = startTime + deltaTime;sv.X=RK*UX;sv.Y=RK*UY;sv.Z=RK*UZ;sv.XDOT=RDOTK*UX+RFDOTK*VX;sv.YDOT=RDOTK*UY+RFDOTK*VY;sv.ZDOT=RDOTK*UZ+RFDOTK*VZ;sv.X = sv.X * XKMPER/AE;sv.Y = sv.Y * XKMPER/AE;sv.Z = sv.Z * XKMPER/AE;sv.XDOT = sv.XDOT * XKMPER/AE*XMNPDA/86400.0;sv.YDOT = sv.YDOT * XKMPER/AE*XMNPDA/86400.0;sv.ZDOT = sv.ZDOT * XKMPER/AE*XMNPDA/86400.0;return sv;}}
/** Copyright (C) 2006 Matthew Funk* Licensed under the Academic Free License version 1.2** 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* Academic Free License version 1.2 for more details.*//*** Implementation of NORAD's SGP4 near-Earth propagator.*/public class SGP4 extends Propagator{private double A1;private double COSIO;private double THETA2;private double X3THM1;private double EOSQ;private double BETAO2;private double BETAO;private double DEL1;private double AO;private double DELO;private double XNODP;private double AODP;private int ISIMP;private double S4;private double QOMS24;private double PERIGE;private double PINVSQ;private double TSI;private double ETA;private double ETASQ;private double EETA;private double PSISQ;private double COEF;private double COEF1;private double c1;private double c2;private double c3;private double c4;private double c5;private double SINIO;private double A3OVK2;private double X1MTH2;private double THETA4;private double TEMP1;private double TEMP2;private double TEMP3;private double XMDOT;private double X1M5TH;private double OMGDOT;private double XHDOT1;private double XNODOT;private double OMGCOF;private double XMCOF;private double XNODCF;private double T2COF;private double XLCOF;private double AYCOF;private double DELMO;private double SINMO;private double X7THM1;private double C1SQ;private double D2;private double D3;private double D4;private double TEMP;private double T3COF;private double T4COF;private double T5COF;private double XMDF;private double OMGADF;private double OMEGA;private double XNODDF;private double XMP;private double TSQ;private double XNODE;private double TEMPA;private double TEMPE;private double TEMPL;private double DELOMG;private double DELM;private double TCUBE;private double TFOUR;private double A;private double E;private double XL;private double BETA;private double XN;private double AXN;private double XLL;private double AYNL;private double XLT;private double AYN;private double CAPU;private double SINEPW;private double COSEPW;private double TEMP4;private double TEMP5;private double TEMP6;private double EPW;private double ECOSE;private double ESINE;private double ELSQ;private double PL;private double R;private double RDOT;private double RFDOT;private double BETAL;private double COSU;private double SINU;private double U;private double COS2U;private double SIN2U;private double RK;private double UK;private double XNODEK;private double XINCK;private double RDOTK;private double RFDOTK;private double SINUK;private double COSUK;private double SINIK;private double COSIK;private double SINNOK;private double COSNOK;private double XMX;private double XMY;private double UX;private double UY;private double UZ;private double VX;private double VY;private double VZ;/*** Constructs an instance of this class.*/public SGP4(String title, ElementSet es) {this.title = title;this.es = es;startTime = es.julian;initialize();}private void initialize() {// RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)// FROM INPUT ELEMENTSA1=Math.pow((XKE/es.XNO),TOTHRD);COSIO=Math.cos(es.XINCL);THETA2=COSIO*COSIO;X3THM1=3.*THETA2-1.;EOSQ=es.EO*es.EO;BETAO2=1.-EOSQ;BETAO=Math.sqrt(BETAO2);DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2);AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)));DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2);XNODP=es.XNO/(1.+DELO);AODP=AO/(1.-DELO);// INITIALIZATION// FOR PERIGEE LESS THAN 220 KILOMETERS, THE ISIMP FLAG IS SET AND// THE EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN SQRT A AND// QUADRATIC VARIATION IN MEAN ANOMALY. ALSO, THE C3 TERM, THE// DELTA OMEGA TERM, AND THE DELTA M TERM ARE DROPPED.ISIMP=0;if ((AODP*(1.-es.EO)/AE) < (220./XKMPER+AE)) ISIMP=1;// FOR PERIGEE BELOW 156 KM, THE VALUES OF// S AND QOMS2T ARE ALTEREDS4=S;QOMS24=QOMS2T;PERIGE=(AODP*(1.-es.EO)-AE)*XKMPER;if (PERIGE < 156.){S4=PERIGE-78.;if (PERIGE <= 98.){S4=20.;}QOMS24=Math.pow(((120.-S4)*AE/XKMPER),4);S4=S4/XKMPER+AE;}PINVSQ=1./(AODP*AODP*BETAO2*BETAO2);TSI=1./(AODP-S4);ETA=AODP*es.EO*TSI;ETASQ=ETA*ETA;EETA=es.EO*ETA;PSISQ=Math.abs(1.-ETASQ);COEF=QOMS24*Math.pow(TSI,4);COEF1=COEF/Math.pow(PSISQ,3.5);c2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)));c1=es.BSTAR*c2;SINIO=Math.sin(es.XINCL);A3OVK2=-XJ3/CK2*Math.pow(AE,3);c3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/es.EO;X1MTH2=1.-THETA2;c4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*(2.+.5*ETASQ)+es.EO*(.5+2.*ETASQ)-2.*CK2*TSI/(AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*(1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*(1.+ETASQ))*Math.cos(2.*es.OMEGAO)));c5=2.*COEF1*AODP*BETAO2*(1.+2.75*(ETASQ+EETA)+EETA*ETASQ);THETA4=THETA2*THETA2;TEMP1=3.*CK2*PINVSQ*XNODP;TEMP2=TEMP1*CK2*PINVSQ;TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP;XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*(13.-78.*THETA2+137.*THETA4);X1M5TH=1.-5.*THETA2;OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4);XHDOT1=-TEMP1*COSIO;XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-7.*THETA2))*COSIO;OMGCOF=es.BSTAR*c3*Math.cos(es.OMEGAO);XMCOF=-TOTHRD*COEF*es.BSTAR*AE/EETA;XNODCF=3.5*BETAO2*XHDOT1*c1;T2COF=1.5*c1;XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO);AYCOF=.25*A3OVK2*SINIO;DELMO=Math.pow((1.+ETA*Math.cos(es.XMO)),3);SINMO=Math.sin(es.XMO);X7THM1=7.*THETA2-1.;if (ISIMP != 1){C1SQ=c1*c1;D2=4.*AODP*TSI*C1SQ;TEMP=D2*TSI*c1/3.;D3=(17.*AODP+S4)*TEMP;D4=.5*TEMP*AODP*TSI*(221.*AODP+31.*S4)*c1;T3COF=D2+2.*C1SQ;T4COF=.25*(3.*D3+c1*(12.*D2+10.*C1SQ));T5COF=.2*(3.*D4+12.*c1*D3+6.*D2*D2+15.*C1SQ*(2.*D2+C1SQ));}}//it appears as though TSINCE is minutespublic OrbitalState propagate(double TSINCE) {// UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAGXMDF=es.XMO+XMDOT*TSINCE;OMGADF=es.OMEGAO+OMGDOT*TSINCE;XNODDF=es.XNODEO+XNODOT*TSINCE;OMEGA=OMGADF;XMP=XMDF;TSQ=TSINCE*TSINCE;XNODE=XNODDF+XNODCF*TSQ;TEMPA=1.-c1*TSINCE;TEMPE=es.BSTAR*c4*TSINCE;TEMPL=T2COF*TSQ;if (ISIMP != 1){DELOMG=OMGCOF*TSINCE;DELM=XMCOF*(Math.pow((1.+ETA*Math.cos(XMDF)),3)-DELMO);TEMP=DELOMG+DELM;XMP=XMDF+TEMP;OMEGA=OMGADF-TEMP;TCUBE=TSQ*TSINCE;TFOUR=TSINCE*TCUBE;TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR;TEMPE=TEMPE+es.BSTAR*c5*(Math.sin(XMP)-SINMO);TEMPL=TEMPL+T3COF*TCUBE+TFOUR*(T4COF+TSINCE*T5COF);}A=AODP*Math.pow(TEMPA,2);E=es.EO-TEMPE;XL=XMP+OMEGA+XNODE+XNODP*TEMPL;BETA=Math.sqrt(1.-E*E);XN=XKE/Math.pow(A,1.5);// LONG PERIOD PERIODICSAXN=E*Math.cos(OMEGA);TEMP=1./(A*BETA*BETA);XLL=TEMP*XLCOF*AXN;AYNL=TEMP*AYCOF;XLT=XL+XLL;AYN=E*Math.sin(OMEGA)+AYNL;// SOLVE KEPLERS EQUATIONCAPU=FMOD2P(XLT-XNODE);TEMP2=CAPU;for (int I = 1 ; I <= 10 ; I++){SINEPW=Math.sin(TEMP2);COSEPW=Math.cos(TEMP2);TEMP3=AXN*SINEPW;TEMP4=AYN*COSEPW;TEMP5=AXN*COSEPW;TEMP6=AYN*SINEPW;EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2;if (Math.abs(EPW-TEMP2) <= E6A) break;TEMP2=EPW;}// SHORT PERIOD PRELIMINARY QUANTITIESECOSE=TEMP5+TEMP6;ESINE=TEMP3-TEMP4;ELSQ=AXN*AXN+AYN*AYN;TEMP=1.-ELSQ;PL=A*TEMP;R=A*(1.-ECOSE);TEMP1=1./R;RDOT=XKE*Math.sqrt(A)*ESINE*TEMP1;RFDOT=XKE*Math.sqrt(PL)*TEMP1;TEMP2=A*TEMP1;BETAL=Math.sqrt(TEMP);TEMP3=1./(1.+BETAL);COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3);SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3);U=ACTAN(SINU,COSU);SIN2U=2.*SINU*COSU;COS2U=2.*COSU*COSU-1.;TEMP=1./PL;TEMP1=CK2*TEMP;TEMP2=TEMP1*TEMP;// UPDATE FOR SHORT PERIODICSRK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U;UK=U-.25*TEMP2*X7THM1*SIN2U;XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U;XINCK=es.XINCL+1.5*TEMP2*COSIO*SINIO*COS2U;RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U;RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1);// ORIENTATION VECTORSSINUK=Math.sin(UK);COSUK=Math.cos(UK);SINIK=Math.sin(XINCK);COSIK=Math.cos(XINCK);SINNOK=Math.sin(XNODEK);COSNOK=Math.cos(XNODEK);XMX=-SINNOK*COSIK;XMY=COSNOK*COSIK;UX=XMX*SINUK+COSNOK*COSUK;UY=XMY*SINUK+SINNOK*COSUK;UZ=SINIK*SINUK;VX=XMX*COSUK-COSNOK*SINUK;VY=XMY*COSUK-SINNOK*SINUK;VZ=SINIK*COSUK;// POSITION AND VELOCITYOrbitalState sv = new OrbitalState();//convert minutes elapsed into fractional days.double deltaTime = TSINCE*60/86400 ;sv.julianTime = startTime + deltaTime - mod;sv.X=RK*UX;sv.Y=RK*UY;sv.Z=RK*UZ;sv.XDOT=RDOTK*UX+RFDOTK*VX;sv.YDOT=RDOTK*UY+RFDOTK*VY;sv.ZDOT=RDOTK*UZ+RFDOTK*VZ;sv.X = sv.X * XKMPER/AE;sv.Y = sv.Y * XKMPER/AE;sv.Z = sv.Z * XKMPER/AE;sv.XDOT = sv.XDOT * XKMPER/AE*XMNPDA/86400.0;sv.YDOT = sv.YDOT * XKMPER/AE*XMNPDA/86400.0;sv.ZDOT = sv.ZDOT * XKMPER/AE*XMNPDA/86400.0;return sv;}}
double toJulian(int year, double epoch) {//use the J2000 reference epochdouble julian = 2451545 + year*365.25 + epoch;return julian;}double getSiderealTime(double JD){double sidereal;double T;double epochTime = JD - 2451545.0;//this is for display -- does not affect calculationsyear = (int)(epochTime / 365.25);double day = epochTime - 365.25*year;year += 2000;int J = (int)(JD-.5);J -=.25;//What?J += 2;int j = J + 32044;int g = j / 146097;int dg = j%146097;int c = (dg/36524 + 1)*3/4;int dc = dg-c*36524;int b = dc/1461;int db = dc % 1461;int a = (db/365 + 1)*3/4;int da = db-a*365;int y = g*400+c*100+b*4+a;int m = (da*5 + 308)/153 - 2;int d = da - (m+4)*153/5+122;Y = y-4800 + (m+2)/12;M = (m+2)% 12 + 1;D = d + 1;double rem = day - (int)day;//half a day for UTC convertif (rem < .5){rem +=.5;day--;}else{rem -= .5;}//fractional leap year? need to look into this. It also occurs in the mod function.if (rem < .25){rem ++;day--;}rem -=.25;//convert to secondsint secs = (int)(rem*86400);hour = secs/3600;secs -= hour*3600;minute = secs/60;secs -= minute*60;second = secs;//end display onlyT = (epochTime) / 36525.0;/* calc mean angle */sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0.000387933 * T * T) - (T * T * T / 38710000.0);/* add a convenient multiple of 360 degrees */sidereal = range_degrees(sidereal);return PI*(sidereal)/180.0;}double range_degrees (double angle){double temp;double range;if (angle >= 0.0 && angle < 360.0)return(angle);temp = (int)(angle / 360);if ( angle < 0.0 )temp --;temp *= 360;range = angle - temp;return (range);}double getCurrentJulianDate(){int JGREG= 15 + 31*(10+12*1582);Date time = new Date();calendar.setTime(time);//get into GMTcalendar.setTimeZone(TimeZone.getTimeZone("GMT"));int year=calendar.get(Calendar.YEAR);int month=calendar.get(Calendar.MONTH) + 1;int day=calendar.get(Calendar.DAY_OF_MONTH);//UT starts at noon, so we need to subtract half a dayhour = calendar.get(Calendar.HOUR_OF_DAY);if (hour < 12){hour += 12;day--;}else{hour -=12;}int julianYear = year;if (year < 0) julianYear++;int julianMonth = month;if (month > 2) {julianMonth++;}else {julianYear--;julianMonth += 13;}double julian = (java.lang.Math.floor(365.25 * julianYear)+ java.lang.Math.floor(30.6001*julianMonth) + day + 1720995.0);if (day + 31 * (month + 12 * year) >= JGREG) {// change over to Gregorian calendarint ja = (int)(0.01 * julianYear);julian += 2 - ja + (0.25 * ja);}//add the fractional timejulian += (hour*3600 + calendar.get(calendar.MINUTE)*60 + calendar.get(calendar.SECOND))/86400.0;return julian;}
class Tuple{float x,y,z;public Tuple(float x, float y, float z){this.x = x;this.y = y;this.z = z;}float dot(Tuple b){return (x*b.x + y*b.y + z*b.z);}void normalize(){float dist = sqrt(x*x + y*y + z*z);x /= dist;y /= dist;z /= dist;}}