diff --git a/default/predict.qth b/default/predict.qth index f393d6f..8963919 100644 --- a/default/predict.qth +++ b/default/predict.qth @@ -2,3 +2,4 @@ W1AW 41.716905 72.727083 25 + 1 diff --git a/predict.c b/predict.c index c9b86d3..6b11389 100644 --- a/predict.c +++ b/predict.c @@ -103,11 +103,11 @@ /* Flow control flag definitions */ #define ALL_FLAGS -1 -#define SGP_INITIALIZED_FLAG 0x000001 /* not used */ +#define SGP_INITIALIZED_FLAG 0x000001 #define SGP4_INITIALIZED_FLAG 0x000002 #define SDP4_INITIALIZED_FLAG 0x000004 -#define SGP8_INITIALIZED_FLAG 0x000008 /* not used */ -#define SDP8_INITIALIZED_FLAG 0x000010 /* not used */ +#define SGP8_INITIALIZED_FLAG 0x000008 +#define SDP8_INITIALIZED_FLAG 0x000010 #define SIMPLE_FLAG 0x000020 #define DEEP_SPACE_EPHEM_FLAG 0x000040 #define LUNAR_TERMS_DONE_FLAG 0x000080 @@ -143,6 +143,7 @@ struct { char callsign[17]; double stnlat; double stnlong; int stnalt; + int mode; } qth; struct { char name[25]; @@ -753,6 +754,148 @@ void select_ephemeris(tle_t *tle) ClearFlag(DEEP_SPACE_EPHEM_FLAG); } +void SGP(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) +{ + /* This function is used to calculate the position and velocity */ + /* of satellites. tsince is time since epoch in minutes, tle is */ + /* a pointer to a tle_t structure with Keplerian orbital */ + /* elements and pos and vel are vector_t structures returning */ + /* ECI satellite position and velocity. Use Convert_Sat_State() */ + /* to convert to km and km/s. */ + + static double ao, qo, xnodot, omgdot, d3o, d2o, d1o, xlo, d4o, c6, c5; + + double c1, c2, c3, c4, cosio, sinio, a1, d1, po, po2no, a, e, p, xnodes, + omgas, xls, axnsl, aynsl, xl, u, eo1, temp5, sineo1, coseo1, + temp2, ecose, esine, el2, pl, pl2, r, rdot, rvdot, temp, sinu, cosu, + su, sin2u, cos2u, rk, uk, xnodek, xinck, sinuk, cosuk, sinnok, + cosnok, sinik, cosik, xmx, xmy, ux, uy, uz, vx, vy, vz; + + int i; + + + /* Initialization */ + + if (isFlagClear(SGP_INITIALIZED_FLAG)) + { + SetFlag(SGP_INITIALIZED_FLAG); + + /* Recover original mean motion (xnodp) and */ + /* semimajor axis (aodp) from input elements. */ + c1=ck2*1.5; + c2=ck2*0.25; + c3=ck2*0.5; + c4=xj3*pow(ae,3)/(4.0*ck2); + cosio=cos(tle->xincl); + sinio=sin(tle->xincl); + a1=pow(xke/tle->xno,tothrd); + d1=c1/a1/a1*(3.*cosio*cosio-1.)/pow(1.-tle->eo*tle->eo,1.5); + ao=a1*(1.-1./3.*d1-d1*d1-134./81.*d1*d1*d1); + po=ao*(1.-tle->eo*tle->eo); + qo=ao*(1.-tle->eo); + xlo=tle->xmo+tle->omegao+tle->xnodeo; + d1o=c3*sinio*sinio; + d2o=c2*(7.*cosio*cosio-1.); + d3o=c1*cosio; + d4o=d3o*sinio; + po2no=tle->xno/(po*po); + omgdot=c1*po2no*(5.*cosio*cosio-1.); + xnodot=-2.*d3o*po2no; + c5=.5*c4*sinio*(3.+5.*cosio)/(1.+cosio); + c6=c4*sinio; + } + + /* Update for secular gravity and atmospheric drag. */ + a=tle->xno+(2.*tle->xndt2o+3.*tle->xndd6o*tsince)*tsince; + a=ao*pow(tle->xno/a,tothrd); + e=e6a; + if (a>qo) + e=1.-qo/a; + p=a*(1.-e*e); + xnodes=tle->xnodeo+xnodot*tsince; + omgas=tle->omegao+omgdot*tsince; + xls=FMod2p(xlo+(tle->xno+omgdot+xnodot+(tle->xndt2o+tle->xndd6o*tsince)*tsince)*tsince); + + /* Long period periodics */ + axnsl=e*cos(omgas); + aynsl=e*sin(omgas)-c6/p; + xl=FMod2p(xls-c5/p*axnsl); + + /* Solve Kepler's Equation */ + u=FMod2p(xl-xnodes); + eo1=u; + temp5=1.; + i=0; + + do + { + sineo1=sin(eo1); + coseo1=cos(eo1); + if (fabs(temp5) < e6a) + break; + temp5=1.-coseo1*axnsl-sineo1*aynsl; + temp5=(u-aynsl*coseo1+axnsl*sineo1-eo1)/temp5; + temp2=fabs(temp5); + if (temp2 > 1.) + temp5=temp2/temp5; + eo1=eo1+temp5; + } while (i++<10); + + /* Short period preliminary quantities */ + ecose=axnsl*coseo1+aynsl*sineo1; + esine=axnsl*sineo1-aynsl*coseo1; + el2=axnsl*axnsl+aynsl*aynsl; + pl=a*(1.-el2); + pl2=pl*pl; + r=a*(1.-ecose); + rdot=xke*sqrt(a)/r*esine; + rvdot=xke*sqrt(pl)/r; + temp=esine/(1.+sqrt(1.-el2)); + sinu=a/r*(sineo1-aynsl-axnsl*temp); + cosu=a/r*(coseo1-axnsl+aynsl*temp); + su=AcTan(sinu,cosu); + + /* Update for short periodics */ + sin2u=(cosu+cosu)*sinu; + cos2u=1.-2.*sinu*sinu; + rk=r+d1o/pl*cos2u; + uk=su-d2o/pl2*sin2u; + xnodek=xnodes+d3o*sin2u/pl2; + xinck=tle->xincl+d4o/pl2*cos2u; + + /* Orientation vectors */ + sinuk=sin(uk); + cosuk=cos(uk); + sinnok=sin(xnodek); + cosnok=cos(xnodek); + sinik=sin(xinck); + cosik=cos(xinck); + 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 velocity */ + pos->x=rk*ux; + pos->y=rk*uy; + pos->z=rk*uz; + vel->x=rvdot*vx+rdot*ux; + vel->y=rvdot*vy+rdot*uy; + vel->z=rvdot*vz+rdot*uz; + + /* Phase in radians */ + phase=xl-xnodes-omgas+twopi; + + if (phase<0.0) + phase+=twopi; + + phase=FMod2p(phase); +} + void SGP4(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) { /* This function is used to calculate the position and velocity */ @@ -1005,6 +1148,302 @@ void SGP4(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) phase=FMod2p(phase); } +void SGP8(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) +{ + /* This function is used to calculate the position and velocity */ + /* of near-earth (period < 225 minutes) satellites. tsince is */ + /* time since epoch in minutes, tle is a pointer to a tle_t */ + /* structure with Keplerian orbital elements and pos and vel */ + /* are vector_t structures returning ECI satellite position and */ + /* velocity. Use Convert_Sat_State() to convert to km and km/s. */ + + static double gamma, ovgpp, ed, qq, xnd, pp, edot, xndt, xnodot, omgdot, + xlldot, xhdt1, xgdt1, xmdot1, unmth2, unm5th, cosio2, sinio2, + sinio, xnodp, tthmun, theta2, cosio; + + + double a1, eosq, betao2, betao, del1, ao, delo, aodp, b, po, pom2, + sing, cosg, temp, theta4, a3cof, pardt1, pardt2, pardt4, tsi, + eta, eta2, psim2, alpha2, eeta, cos2g, d5, d1, d2, d3, d4, + b1, b2, b3, co, c1, c4, c5, xndtn, d6, d7, d8, c8, c9, d2o, aldtal, + tsdtts, etdt, psdtps, sin2g, codtco, c1dtc1, d9, d1o, d11, d12, + d13, d14, d15, d1dt, d2dt, d3dt, d4dt, d5dt, c4dt, c5dt, d16, + xnddt, eddot, d25, d17, tsddts, etddt, d18, d19, d23, d1ddt, + xntrdt, tmnddt, xmam, omgasm, xnodes, xn, em, z1, temp1, z7, zc2, + sine, cose, zc5, cape, am, beta2m, sinos, cosos, axnm, aynm, pm, + g1, g2, g3, beta, g4, g5, snf, csf, snfg, csfg, sn2f2g, cs2f2g, + ecosf, g1o, rm, aovr, g13, g14, dr, diwc, di, sin2du, xlamb, y4, + y5, r, rdot, rvdot, snlamb, cslamb, ux, vx, uy, vy, uz, vz; + + int i; + + /* Initialization */ + + if (isFlagClear(SGP8_INITIALIZED_FLAG)) + { + SetFlag(SGP8_INITIALIZED_FLAG); + + /* Recover original mean motion (xnodp) and */ + /* semimajor axis (aodp) from input elements. */ + + a1=pow(xke/tle->xno,tothrd); + cosio=cos(tle->xincl); + theta2=cosio*cosio; + tthmun=3.*theta2-1.; + eosq=tle->eo*tle->eo; + betao2=1.-eosq; + betao=sqrt(betao2); + del1=1.5*ck2*tthmun/(a1*a1*betao*betao2); + ao=a1*(1.-del1*(.5*tothrd+del1*(1.+134./81.*del1))); + delo=1.5*ck2*tthmun/(ao*ao*betao*betao2); + aodp=ao/(1.-delo); + xnodp=tle->xno/(1.+delo); + b=2.*tle->bstar/rho; + + /* Initialization */ + ClearFlag(SIMPLE_FLAG); + po=aodp*betao2; + pom2=1./(po*po); + sinio=sin(tle->xincl); + sing=sin(tle->omegao); + cosg=cos(tle->omegao); + temp=.5*tle->xincl; + sinio2=sin(temp); + cosio2=cos(temp); + theta4=pow(theta2,2); + unm5th=1.-5.*theta2; + unmth2=1.-theta2; + a3cof=-xj3/ck2*pow(ae,3); + pardt1=3.*ck2*pom2*xnodp; + pardt2=pardt1*ck2*pom2; + pardt4=1.25*ck4*pom2*pom2*xnodp; + xmdot1=.5*pardt1*betao*tthmun; + xgdt1=-.5*pardt1*unm5th; + xhdt1=-pardt1*cosio; + xlldot=xnodp+xmdot1+.0625*pardt2*betao*(13.-78.*theta2+137.*theta4); + omgdot=xgdt1+.0625*pardt2*(7.-114.*theta2+395.*theta4)+pardt4*(3.-36.*theta2+49.*theta4); + xnodot=xhdt1+(.5*pardt2*(4.-19.*theta2)+2.*pardt4*(3.-7.*theta2))*cosio; + tsi=1./(po-s); + eta=tle->eo*s*tsi; + eta2=pow(eta,2); + psim2=fabs(1./(1.-eta2)); + alpha2=1.+eosq; + eeta=tle->eo*eta; + cos2g=2.*cosg*cosg-1; + d5=tsi*psim2; + d1=d5/po; + d2=12.+eta2*(36.+4.5*eta2); + d3=eta2*(15.+2.5*eta2); + d4=eta*(5.+3.75*eta2); + b1=ck2*tthmun; + b2=-ck2*unmth2; + b3=a3cof*sinio; + co=.5*b*rho*qoms2t*xnodp*aodp*pow(tsi,4)*pow(psim2,3.5)/sqrt(alpha2); + c1=1.5*xnodp*pow(alpha2,2)*co; + c4=d1*d3*b2; + c5=d5*d4*b3; + xndt=c1*((2.+eta2*(3.+34.*eosq)+5.*eeta*(4.+eta2)+8.5*eosq)+d1*d2*b1+c4*cos2g+c5*sing); + xndtn=xndt/xnodp; + + /* If drag is very small, the SIMPLE_FLAG is set and */ + /* the equations are truncated to linear variation */ + /* in mean motion and quadratic variation in mean */ + /* anomaly. */ + if (fabs(xndt*xmnpda)<2.16E-3) + { + SetFlag(SIMPLE_FLAG); + edot=-tothrd*xndtn*(1.-tle->eo); + } + else + { + d6=eta*(30.+22.5*eta2); + d7=eta*(5.+12.5*eta2); + d8=1.+eta2*(6.75+eta2); + c8=d1*d7*b2; + c9=d5*d8*b3; + edot=-co*(eta*(4.+eta2+eosq*(15.5+7.*eta2))+tle->eo*(5.+15.*eta2)+d1*d6*b1+c8*cos2g+c9*sing); + d2o=.5*tothrd*xndtn; + aldtal=tle->eo*edot/alpha2; + tsdtts=2.*aodp*tsi*(d2o*betao2+tle->eo*edot); + etdt=(edot+tle->eo*tsdtts)*tsi*s; + psdtps=-eta*etdt*psim2; + sin2g=2.*sing*cosg; + codtco=d2o+4.*tsdtts-aldtal-7.*psdtps; + c1dtc1=xndtn+4.*aldtal+codtco; + d9=eta*(6.+68.*eosq)+tle->eo*(20.+15.*eta2); + d1o=.5*eta*(4.+eta2)+tle->eo*(17.+68.*eta2); + d11=eta*(72.+18.*eta2); + d12=eta*(30.+10.*eta2); + d13=5.+11.25*eta2; + d14=tsdtts-2.*psdtps; + d15=2.*(d2o+tle->eo*edot/betao2); + d1dt=d1*(d14+d15); + d2dt=etdt*d11; + d3dt=etdt*d12; + d4dt=etdt*d13; + d5dt=d5*d14; + c4dt=b2*(d1dt*d3+d1*d3dt); + c5dt=b3*(d5dt*d4+d5*d4dt); + d16= + d9*etdt+d1o*edot+ + b1*(d1dt*d2+d1*d2dt)+ + c4dt*cos2g+c5dt*sing+xgdt1*(c5*cosg-2.*c4*sin2g); + xnddt=c1dtc1*xndt+c1*d16; + eddot=codtco*edot-co*( + (4.+3.*eta2+30.*eeta+eosq*(15.5+21.*eta2))*etdt+(5.+15.*eta2 + +eeta*(31.+14.*eta2))*edot+ + b1*(d1dt*d6+d1*etdt*(30.+67.5*eta2))+ + b2*(d1dt*d7+d1*etdt*(5.+37.5*eta2))*cos2g+ + b3*(d5dt*d8+d5*etdt*eta*(13.5+4.*eta2))*sing+xgdt1*(c9* + cosg-2.*c8*sin2g)); + d25=edot*edot; + d17=xnddt/xnodp-xndtn*xndtn; + tsddts=2.*tsdtts*(tsdtts-d2o)+aodp*tsi*(tothrd*betao2*d17-4.*d2o* + tle->eo*edot+2.*(d25+tle->eo*eddot)); + etddt=(eddot+2.*edot*tsdtts)*tsi*s+tsddts*eta; + d18=tsddts-tsdtts*tsdtts; + d19=-psdtps*psdtps/eta2-eta*etddt*psim2-psdtps*psdtps; + d23=etdt*etdt; + d1ddt=d1dt*(d14+d15)+d1*(d18-2.*d19+tothrd*d17+2.*(alpha2*d25 + /betao2+tle->eo*eddot)/betao2); + xntrdt=xndt*(2.*tothrd*d17+3.* + (d25+tle->eo*eddot)/alpha2-6.*aldtal*aldtal+ + 4.*d18-7.*d19)+ + c1dtc1*xnddt+c1*(c1dtc1*d16+ + d9*etddt+d1o*eddot+d23*(6.+30.*eeta+68.*eosq)+ + etdt*edot*(40.+30.* + eta2+272.*eeta)+d25*(17.+68.*eta2)+ + b1*(d1ddt*d2+2.*d1dt*d2dt+d1*(etddt*d11+d23*(72.+54.*eta2)))+ + b2*(d1ddt*d3+2.*d1dt*d3dt+d1*(etddt*d12+d23*(30.+30.*eta2)))* + cos2g+ + b3*((d5dt*d14+d5*(d18-2.*d19))* + d4+2.*d4dt*d5dt+d5*(etddt*d13+22.5*eta*d23))*sing+xgdt1* + ((7.*d2o+4.*tle->eo*edot/betao2)* + (c5*cosg-2.*c4*sin2g) + +((2.*c5dt*cosg-4.*c4dt*sin2g)-xgdt1*(c5*sing+4.* + c4*cos2g)))); + tmnddt=xnddt*1.E9; + temp=tmnddt*tmnddt-xndt*1.E18*xntrdt; + pp=(temp+tmnddt*tmnddt)/temp; + gamma=-xntrdt/(xnddt*(pp-2.)); + xnd=xndt/(pp*gamma); + qq=1.-eddot/(edot*gamma); + ed=edot/(qq*gamma); + ovgpp=1./(gamma*(pp+1.)); + } + } + + /* Update for secular gravity and atmospheric drag */ + xmam=FMod2p(tle->xmo+xlldot*tsince); + omgasm=tle->omegao+omgdot*tsince; + xnodes=tle->xnodeo+xnodot*tsince; + if (isFlagSet(SIMPLE_FLAG)) + { + xn=xnodp+xndt*tsince; + em=tle->eo+edot*tsince; + z1=.5*xndt*tsince*tsince; + } + else + { + temp=1.-gamma*tsince; + temp1=pow(temp,pp); + xn=xnodp+xnd*(1.-temp1); + em=tle->eo+ed*(1.-pow(temp,qq)); + z1=xnd*(tsince+ovgpp*(temp*temp1-1.)); + } + z7=3.5*tothrd*z1/xnodp; + xmam=FMod2p(xmam+z1+z7*xmdot1); + omgasm=omgasm+z7*xgdt1; + xnodes=xnodes+z7*xhdt1; + + /* Solve Kepler's Equation */ + zc2=xmam+em*sin(xmam)*(1.+em*cos(xmam)); + i=0; + + do + { + sine=sin(zc2); + cose=cos(zc2); + zc5=1./(1.-em*cose); + cape=(xmam+em*sine-zc2)*zc5+zc2; + if (fabs(cape-zc2)<=e6a) break; + zc2=cape; + } while (i++<10); + + /* Short period preliminary quantities */ + am=pow(xke/xn,tothrd); + beta2m=1.-em*em; + sinos=sin(omgasm); + cosos=cos(omgasm); + axnm=em*cosos; + aynm=em*sinos; + pm=am*beta2m; + g1=1./pm; + g2=.5*ck2*g1; + g3=g2*g1; + beta=sqrt(beta2m); + g4=.25*a3cof*sinio; + g5=.25*a3cof*g1; + snf=beta*sine*zc5; + csf=(cose-em)*zc5; + fm=AcTan(snf,csf); + snfg=snf*cosos+csf*sinos; + csfg=csf*cosos-snf*sinos; + sn2f2g=2.*snfg*csfg; + cs2f2g=2.*csfg*csfg-1.; + ecosf=em*csf; + g1o=fm-xmam+em*snf; + rm=pm/(1.+ecosf); + aovr=am/rm; + g13=xn*aovr; + g14=-g13*aovr; + dr=g2*(unmth2*cs2f2g-3.*tthmun)-g4*snfg; + diwc=3.*g3*sinio*cs2f2g-g5*aynm; + di=diwc*cosio; + + /* Update for short period periodics */ + sin2du=sinio2*( + g3*(.5*(1.-7.*theta2)*sn2f2g-3.*unm5th*g1o)-g5*sinio*csfg*(2.+ + ecosf))-.5*g5*theta2*axnm/cosio2; + xlamb=fm+omgasm+xnodes+g3*(.5*(1.+6.*cosio-7.*theta2)*sn2f2g-3.* + (unm5th+2.*cosio)*g1o)+g5*sinio*(cosio*axnm/(1.+cosio)-(2. + +ecosf)*csfg); + y4=sinio2*snfg+csfg*sin2du+.5*snfg*cosio2*di; + y5=sinio2*csfg-snfg*sin2du+.5*csfg*cosio2*di; + r=rm+dr; + rdot=xn*am*em*snf/beta+g14*(2.*g2*unmth2*sn2f2g+g4*csfg); + rvdot=xn*am*am*beta/rm+g14*dr+am*g13*sinio*diwc; + + /* Orientation vectors */ + snlamb=sin(xlamb); + cslamb=cos(xlamb); + temp=2.*(y5*snlamb-y4*cslamb); + ux=y4*temp+cslamb; + vx=y5*temp-snlamb; + temp=2.*(y5*cslamb+y4*snlamb); + uy=-y4*temp+snlamb; + vy=-y5*temp+cslamb; + temp=2.*sqrt(1.-y4*y4-y5*y5); + uz=y4*temp; + vz=y5*temp; + + /* Position and velocity */ + pos->x=r*ux; + pos->y=r*uy; + pos->z=r*uz; + vel->x=rdot*ux+rvdot*vx; + vel->y=rdot*uy+rvdot*vy; + vel->z=rdot*uz+rvdot*vz; + + /* Phase in radians */ + phase=xlamb-xnodes-omgasm+twopi; + + if (phase<0.0) + phase+=twopi; + + phase=FMod2p(phase); +} + void Deep(int ientry, tle_t * tle, deep_arg_t * deep_arg) { /* This function is used by SDP4 to add lunar and solar */ @@ -1750,6 +2189,218 @@ void SDP4(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) phase=FMod2p(phase); } +void SDP8(double tsince, tle_t * tle, vector_t * pos, vector_t * vel) +{ + /* This function is used to calculate the position and velocity */ + /* of deep-space (period > 225 minutes) satellites. tsince is */ + /* time since epoch in minutes, tle is a pointer to a tle_t */ + /* structure with Keplerian orbital elements and pos and vel */ + /* are vector_t structures returning ECI satellite position and */ + /* velocity. Use Convert_Sat_State() to convert to km and km/s. */ + + int i; + + static double edot, xndot, xmdot1, xhdot1, xgdot1, unmth2, unm5th, + cosio2, sinio2, tthmun; + + double a1, del1, ao, delo, b, po, pom2, temp, theta4, a3cof, pardt1, + pardt2, pardt4, tsi, eta, eta2, psim2, alpha2, eeta, cos2g, + d5, d1, d2, d3, d4, b1, b2, b3, co, c1, c4, c5, xndotn, z1, + z7, zc2, sine, cose, zc5, cape, am, beta2m, sinos, cosos, + axnm, aynm, pm, g1, g2, g3, beta, g4, g5, snf, csf, snfg, + csfg, sn2f2g, cs2f2g, ecosf, g1o, rm, aovr, g13, g14, dr, + diwc, di, sini2, sni2du, xlamb, y4, y5, r, rdot, rvdot, + snlamb, cslamb, ux, vx, uy, vy, uz, vz; + + static deep_arg_t deep_arg; + + /* Initialization */ + + if (isFlagClear(SDP8_INITIALIZED_FLAG)) + { + SetFlag(SDP8_INITIALIZED_FLAG); + + /* Recover original mean motion (xnodp) and */ + /* semimajor axis (aodp) from input elements. */ + /* Calculate ballistic coefficient (b term) */ + /* from input b* drag term */ + a1=pow(xke/tle->xno,tothrd); + deep_arg.cosio=cos(tle->xincl); + deep_arg.theta2=deep_arg.cosio*deep_arg.cosio; + tthmun=3.*deep_arg.theta2-1.; + deep_arg.eosq=tle->eo*tle->eo; + deep_arg.betao2=1.-deep_arg.eosq; + deep_arg.betao=sqrt(deep_arg.betao2); + del1=1.5*ck2*tthmun/(a1*a1*deep_arg.betao*deep_arg.betao2); + ao=a1*(1.-del1*(.5*tothrd+del1*(1.+134./81.*del1))); + delo=1.5*ck2*tthmun/(ao*ao*deep_arg.betao*deep_arg.betao2); + deep_arg.aodp=ao/(1.-delo); + deep_arg.xnodp=tle->xno/(1.+delo); + b=2.*tle->bstar/rho; + + /* Initialization */ + po=deep_arg.aodp*deep_arg.betao2; + pom2=1./(po*po); + deep_arg.sinio=sin(tle->xincl); + deep_arg.sing=sin(tle->omegao); + deep_arg.cosg=cos(tle->omegao); + temp=.5*tle->xincl; + sinio2=sin(temp); + cosio2=cos(temp); + theta4=deep_arg.theta2*deep_arg.theta2; + unm5th=1.-5.*deep_arg.theta2; + unmth2=1.-deep_arg.theta2; + a3cof=-xj3/ck2*ae*ae*ae; + pardt1=3.*ck2*pom2*deep_arg.xnodp; + pardt2=pardt1*ck2*pom2; + pardt4=1.25*ck4*pom2*pom2*deep_arg.xnodp; + xmdot1=.5*pardt1*deep_arg.betao*tthmun; + xgdot1=-.5*pardt1*unm5th; + xhdot1=-pardt1*deep_arg.cosio; + deep_arg.xmdot=deep_arg.xnodp+xmdot1+ + .0625*pardt2*deep_arg.betao*(13.-78.*deep_arg.theta2+137.*theta4); + deep_arg.omgdot=xgdot1+ + .0625*pardt2*(7.-114.*deep_arg.theta2+395.*theta4)+pardt4*(3.-36.* + deep_arg.theta2+49.*theta4); + deep_arg.xnodot=xhdot1+ + (.5*pardt2*(4.-19.*deep_arg.theta2)+2.*pardt4*(3.-7.*deep_arg.theta2))*deep_arg.cosio; + tsi=1./(po-s); + eta=tle->eo*s*tsi; + eta2=eta*eta; + psim2=fabs(1./(1.-eta2)); + alpha2=1.+deep_arg.eosq; + eeta=tle->eo*eta; + cos2g=2.*deep_arg.cosg*deep_arg.cosg-1.; + d5=tsi*psim2; + d1=d5/po; + d2=12.+eta2*(36.+4.5*eta2); + d3=eta2*(15.+2.5*eta2); + d4=eta*(5.+3.75*eta2); + b1=ck2*tthmun; + b2=-ck2*unmth2; + b3=a3cof*deep_arg.sinio; + co=.5*b*rho*qoms2t*deep_arg.xnodp*deep_arg.aodp*pow(tsi,4)*pow(psim2,3.5)/sqrt(alpha2); + c1=1.5*deep_arg.xnodp*alpha2*alpha2*co; + c4=d1*d3*b2; + c5=d5*d4*b3; + xndot=c1*( + (2.+eta2*(3.+34.*deep_arg.eosq)+5.*eeta*(4.+eta2)+8.5*deep_arg.eosq)+ + d1*d2*b1+c4*cos2g+c5*deep_arg.sing); + xndotn=xndot/deep_arg.xnodp; + edot=-tothrd*xndotn*(1.-tle->eo); + + /* initialize Deep() */ + + Deep(dpinit,tle,&deep_arg); + } + + /* Update for secular gravity and atmospheric drag */ + z1=.5*xndot*tsince*tsince; + z7=3.5*tothrd*z1/deep_arg.xnodp; + deep_arg.xll=tle->xmo+deep_arg.xmdot*tsince; + deep_arg.omgadf=tle->omegao+deep_arg.omgdot*tsince+z7*xgdot1; + deep_arg.xnode=tle->xnodeo+deep_arg.xnodot*tsince+z7*xhdot1; + deep_arg.xn=deep_arg.xnodp; + deep_arg.t=tsince; + Deep(dpsec,tle,&deep_arg); + deep_arg.xn=deep_arg.xn+xndot*tsince; + deep_arg.em=deep_arg.em+edot*tsince; + deep_arg.xll=deep_arg.xll+z1+z7*xmdot1; + Deep(dpper,tle,&deep_arg); + deep_arg.xll=FMod2p(deep_arg.xll); + + /* Solve Kepler's Equation */ + zc2=deep_arg.xll+deep_arg.em*sin(deep_arg.xll)*(1.+deep_arg.em*cos(deep_arg.xll)); + i=0; + + do + { + sine=sin(zc2); + cose=cos(zc2); + zc5=1./(1.-deep_arg.em*cose); + cape=(deep_arg.xll+deep_arg.em*sine-zc2)* + zc5+zc2; + if (fabs(cape-zc2)<=e6a) + break; + zc2=cape; + } while (i++<10); + + /* Short period preliminary quantities */ + am=pow(xke/deep_arg.xn,tothrd); + beta2m=1.-deep_arg.em*deep_arg.em; + sinos=sin(deep_arg.omgadf); + cosos=cos(deep_arg.omgadf); + axnm=deep_arg.em*cosos; + aynm=deep_arg.em*sinos; + pm=am*beta2m; + g1=1./pm; + g2=.5*ck2*g1; + g3=g2*g1; + beta=sqrt(beta2m); + g4=.25*a3cof*deep_arg.sinio; + g5=.25*a3cof*g1; + snf=beta*sine*zc5; + csf=(cose-deep_arg.em)*zc5; + fm=AcTan(snf,csf); + snfg=snf*cosos+csf*sinos; + csfg=csf*cosos-snf*sinos; + sn2f2g=2.*snfg*csfg; + cs2f2g=2.*csfg*csfg-1.; + ecosf=deep_arg.em*csf; + g1o=fm-deep_arg.xll+deep_arg.em*snf; + rm=pm/(1.+ecosf); + aovr=am/rm; + g13=deep_arg.xn*aovr; + g14=-g13*aovr; + dr=g2*(unmth2*cs2f2g-3.*tthmun)-g4*snfg; + diwc=3.*g3*deep_arg.sinio*cs2f2g-g5*aynm; + di=diwc*deep_arg.cosio; + sini2=sin(.5*deep_arg.xinc); + + /* Update for short periodics */ + sni2du=sinio2*( + g3*(.5*(1.-7.*deep_arg.theta2)*sn2f2g-3.*unm5th*g1o)-g5*deep_arg.sinio*csfg*(2.+ + ecosf))-.5*g5*deep_arg.theta2*axnm/cosio2; + xlamb=fm+deep_arg.omgadf+deep_arg.xnode+g3*(.5*(1.+6.*deep_arg.cosio-7.*deep_arg.theta2)*sn2f2g-3.* + (unm5th+2.*deep_arg.cosio)*g1o)+g5*deep_arg.sinio*(deep_arg.cosio*axnm/(1.+deep_arg.cosio)-(2. + +ecosf)*csfg); + y4=sini2*snfg+csfg*sni2du+.5*snfg*cosio2*di; + y5=sini2*csfg-snfg*sni2du+.5*csfg*cosio2*di; + r=rm+dr; + rdot=deep_arg.xn*am*deep_arg.em*snf/beta+g14*(2.*g2*unmth2*sn2f2g+g4*csfg); + rvdot=deep_arg.xn*am*am*beta/rm+ + g14*dr+am*g13*deep_arg.sinio*diwc; + + /* Orientation vectors */ + snlamb=sin(xlamb); + cslamb=cos(xlamb); + temp=2.*(y5*snlamb-y4*cslamb); + ux=y4*temp+cslamb; + vx=y5*temp-snlamb; + temp=2.*(y5*cslamb+y4*snlamb); + uy=-y4*temp+snlamb; + vy=-y5*temp+cslamb; + temp=2.*sqrt(1.-y4*y4-y5*y5); + uz=y4*temp; + vz=y5*temp; + + /* Position and velocity */ + pos->x=r*ux; + pos->y=r*uy; + pos->z=r*uz; + vel->x=rdot*ux+rvdot*vx; + vel->y=rdot*uy+rvdot*vy; + vel->z=rdot*uz+rvdot*vz; + + /* Phase in radians */ + phase=xlamb-deep_arg.xnode-deep_arg.omgadf+twopi; + + if (phase<0.0) + phase+=twopi; + + phase=FMod2p(phase); +} + void Calculate_User_PosVel(double time, geodetic_t *geodetic, vector_t *obs_pos, vector_t *obs_vel) { /* Calculate_User_PosVel() passes the user's geodetic position @@ -2719,6 +3370,7 @@ char ReadDataFiles() fscanf(fd,"%lf", &qth.stnlat); fscanf(fd,"%lf", &qth.stnlong); fscanf(fd,"%d", &qth.stnalt); + fscanf(fd,"%d", &qth.mode); fclose(fd); obs_geodetic.lat=qth.stnlat*deg2rad; @@ -2953,6 +3605,10 @@ void SaveQTH() { /* This function saves QTH data to the QTH data file. */ + //Catch a bad mode + if(qth.mode>2||qth.mode<0) + qth.mode = 1; + FILE *fd; fd=fopen(qthfile,"w"); @@ -2961,6 +3617,7 @@ void SaveQTH() fprintf(fd," %g\n",qth.stnlat); fprintf(fd," %g\n",qth.stnlong); fprintf(fd," %d\n",qth.stnalt); + fprintf(fd," %d\n",qth.mode); fclose(fd); } @@ -3761,18 +4418,40 @@ void Calc() age=jul_utc-jul_epoch; /* Copy the ephemeris type in use to ephem string. */ + /* Call NORAD routines according to deep-space flag. */ - if (isFlagSet(DEEP_SPACE_EPHEM_FLAG)) + if (qth.mode==0) //SGP + { + strcpy(ephem,"SGP"); + SGP(tsince, &tle, &pos, &vel); + } + else if (qth.mode==1) //SxP4 + { + if(isFlagSet(DEEP_SPACE_EPHEM_FLAG)) + { strcpy(ephem,"SDP4"); + SDP4(tsince, &tle, &pos, &vel); + } else + { strcpy(ephem,"SGP4"); + SGP4(tsince, &tle, &pos, &vel); + } + } + else if (qth.mode==2) //SxP8 + { + if(isFlagSet(DEEP_SPACE_EPHEM_FLAG)) + { + strcpy(ephem,"SDP8"); + SDP8(tsince, &tle, &pos, &vel); + } + else + { + strcpy(ephem,"SGP8"); + SGP8(tsince, &tle, &pos, &vel); + } + } - /* Call NORAD routines according to deep-space flag. */ - - if (isFlagSet(DEEP_SPACE_EPHEM_FLAG)) - SDP4(tsince, &tle, &pos, &vel); - else - SGP4(tsince, &tle, &pos, &vel); /* Scale position and velocity vectors to km and km/sec */ @@ -4791,25 +5470,28 @@ void QthEdit() mvprintw(7,0,"\t\t * Ground Station Location Editing Utility *\n\n\n"); attrset(COLOR_PAIR(4)|A_BOLD); - printw("\n\t\t\tStation Callsign :"); - printw("\n\t\t\tStation Latitude :"); - printw("\n\t\t\tStation Longitude :"); - printw("\n\t\t\tStation Altitude :"); + printw("\n\t\t\tStation Callsign :"); + printw("\n\t\t\tStation Latitude :"); + printw("\n\t\t\tStation Longitude :"); + printw("\n\t\t\tStation Altitude :"); + printw("\n\t\t\tPerturbation Model :"); attrset(COLOR_PAIR(2)|A_BOLD); - mvprintw(11,44,"%s",qth.callsign); + mvprintw(11,45,"%s",qth.callsign); if (io_lat=='N') - mvprintw(12,44,"%g [DegN]",+qth.stnlat); + mvprintw(12,45,"%g [DegN]",+qth.stnlat); else - mvprintw(12,44,"%g [DegS]",-qth.stnlat); + mvprintw(12,45,"%g [DegS]",-qth.stnlat); if (io_lon=='W') - mvprintw(13,44,"%g [DegW]",+qth.stnlong); + mvprintw(13,45,"%g [DegW]",+qth.stnlong); else - mvprintw(13,44,"%g [DegE]",-qth.stnlong); + mvprintw(13,45,"%g [DegE]",-qth.stnlong); - mvprintw(14,44,"%d [m]",qth.stnalt); + mvprintw(14,45,"%d [m]",qth.stnalt); + + mvprintw(15,45,"%d",qth.mode); refresh(); @@ -4817,7 +5499,7 @@ void QthEdit() mvprintw(18,12,"Enter the callsign or identifier of your ground station"); - if (KbEdit(45,12)) + if (KbEdit(46,12)) strncpy(qth.callsign,temp,16); if (io_lat=='N') @@ -4832,7 +5514,7 @@ void QthEdit() mvprintw(19,12," Decimal (74.2467) or DMS (74 14 48) format allowed"); - if (KbEdit(45,13)) + if (KbEdit(46,13)) { if (io_lat=='N') qth.stnlat=+ReadBearing(temp); @@ -4850,7 +5532,7 @@ void QthEdit() else mvprintw(18,12,"Enter your longitude in degrees EAST (west=negative) "); - if (KbEdit(45,14)) + if (KbEdit(46,14)) { if (io_lon=='W') qth.stnlong=+ReadBearing(temp); @@ -4864,9 +5546,20 @@ void QthEdit() sprintf(temp,"%d",qth.stnalt); - if (KbEdit(45,15)) + if (KbEdit(46,15)) sscanf(temp,"%d",&qth.stnalt); + move(18,12); + clrtoeol(); + + sprintf(temp,"%d",qth.mode); + + mvprintw(18,20, " Enter the perturbation model to use."); + mvprintw(19,20, "0 = SGP 1 = SGP4/SDP4 2 = SGP8/SDP8"); + + if (KbEdit(46,16)) + sscanf(temp, "%d",&qth.mode); + if (resave) { SaveQTH();