This finds RA DEC and distance earth to comet given 6 orbit elements an X,Y,Z sun, ie, and k Elements dt = 40 days past perihelion for this example eccentricity = ee = 0.2 perihelion = qq = 0.4244 AU inclination= ii = 72 deg argument of perihelion= ww = 105 deg Longitude of Ascending node = om = 293 deg inclination of earth = ie = 23.441028 deg Planetary gravity contant = kk = 0.01720209895 Sun position for this example, earth center system xs = -0.931108260968 AU ys = 0.371439715781 AU zs = 0.161052202235 AU ========= sqrt is square root symbol |> is right triangle symbol theta is theta symbol radec1() :Prgm :ClIO :40->dt days :xs = -0.931108260968 AU :ys = 0.371439715781 AU :zs = 0.161052202235 AU :0.01720209895 -> kk :0.2 -> ee :0.4255 -> qq AU :72 -> ii degrees :105 -> ww degrees :293 -> om degrees :23.441028 -> ie degrees :SetMode("Angle","Radian") :qq/(1-ee) -> aa Semimajor axis = 0.531875 AU :kk*dt/aa^(1.5) -> mm Mean Anomaly = 1.77389155706 radians :zeros(ea-ee*sin(e1)-mm,e1=mm) -> dddy Kepler's Equation :Product(dddy) -> e1 e1 = eccentric anomaly 1.95900897925 radians :aa*(cos(e1)-ee) -> x x = -0.307708130153 AU :aa* sqrt(1-ee^2)*sin(e1) -> y y = 0.492350232585 AU :Setmode("Angle","Degree") :R|>Ptheta(x,y) -> th th = true anomaly = 122.535231561 deg :R|>Pr(x,y) -> rr radius = 0.572141626031 AU :P|>Ex(rr,th+ww) -> xx1 xx1 = -0.386273822554 rotate around ww :P|>Ry(rr,th+ww) -> yy0 yy0 = -0.422064656477 :P|>Rx(yy0,ii) -> yy1 yy1 = -0.130425151576 rotate around ii :P|>Ry(yy0,ii) -> zz1 zz = -0.40140734184 :R|>Ptheta(xx1,yy1) -> tt tt = -161.342773319 :R|>Pr(xx1,yy1) -> rom rom = 131.657226681 :P|>Rx(rom,tt+om) -> xx3 xx3 = -0.270986191631 rotate around om :P|>Ry(rom,tt+om) -> yy2 yy2 = 0.304605761763 :R|>Ptheta(yy2,zz1)->tha tha = -52.8071544636 :R|>Pr(yy2,zz1) -> ram ram = 0.503897334963 :P|>Rx(ram,(tha+ie) -> yy3 yy3 = 0.439148484293 rotate around ie :xx3+xs->ss1 ss1 = -1.2020944526 :yy3+ys->ss2 ss2 = 0.810588200074 :zz3+zs -> ss3 ss3= -0.086053307469 :sqrt(ss1*ss1+ss2*ss2+ss3*ss3) -> r0 r0= 1.45240816398 Earth to comet AU :Disp "r0",r0 :sin-1(ss3/rr0) -> dc :Disp "dc",dc dc = -3.39669019616 Declination degrees :R|>Ptheta(ss1,ss2) -> ra ra= 146.007690781 degrees RA :Disp "ra",ra :EndProgm ======== Define true elements ra0() Prgm SetMode("Exact/Approx","Approximate") 0.2->ee inclination of orbit degrees 0.4244->qq perihelion in AU 72->ii inclination of orbit degrees 105->ww argument of perigee degrees 293->om longitude of the ascending node degrees 23.441028->ie tilt of earth degrees 0.01720209895->kk gravity constant EndPrgm theta is theta symbol sqrt is square root symbol values are for time dt=40 days Needs kk,qq,ee,ii,om,w,ie,dt,xs,ys,zs :radec() :Prgm :ClIO :SetMode(Angle","Radian") :SetMode("Exact/Approx","Approximate") :qq/(1-ee) -> aa Semimajor axis = 0.531875 AU :kk*dt/aa^(1.5) -> mm Mean Anomaly = 1.77389155706 readians :zeros(ea-ee*sin(e1)-mm,e1) -> dddy :Product(dddy) -> e1 e1 = eccentric anomaly 1.95900897925 radians :aa * (cos(e1)-ee) -> x x = -0.307708130153 :aa * sqrt(1-ee^2)*sin(e1) -> y y = 0.492350232585 :Setmode("Angle","Degree") :R|>Ptheta(x,y) -> th th = true anomaly = 122.535231561 deg :R|>Pr(x,y) -> rr radius = 0.572141626031 AU :P|>Ex(rr,th+ww) -> xx1 xx1 = -0.386273822554 rotate around ww :P|>Ry(rr,th+ww) -> yy0 yy0 = -0.422064656477 AU :P|>Rx(yy0,ii) -> yy1 yy1 = -0.130425151576 rotate around ii :P|>Ry(yy0,ii) -> zz zz = -0.40140734184 AU :R|>Ptheta(xx1,yy1) -> tt tt = -161.342773319 degrees :R|>Pr(xx1,yy1) -> rom rom = 0.40769864625 AU :P|>Rx(rom,tt+om) -> xx3 xx3 = -0.270986191631 AU rotate around om :P|>Ry(rom,tt+om) -> yy2 yy2 = 0.304605761763 AU :R|>Ptheta(yy2,zz)->tha tha = -52.8071544636 degrees :R|>Pr(yy2,zz) -> ram ram = 0.503897334963 AU :P|>Rx(ram,(tha+ie) -> yy3 yy3 = 0.439148484293 AU rotate around ie :xx3+xs -> ss1 ss1 = -1.2020944526 AU :yy3+ys -> ss2 ss2 = 0.810588200074 AU :zz3+zs -> ss3 ss3= -0.086053307469 AU :sqrt(ss1*ss1+ss2*ss2+ss3*ss3) -> r0 r0= 1.45240816398 AU Earth to comet AU :Disp "r0",r0 :sin-1(ss3/rr0) -> dc :Disp "dc",dc dc = -3.39669019616 Declination degrees :R|>Ptheta(ss1,ss2) -> ra ra= 146.007690781 Right Ascension degree :Disp "ra",ra :EndProgm ----------- :ra11() time 1=35 days :Prgm :ClIO :35->dt :dt->dt1 :-0.893557086283-xs :0.441098630334->ys :0.191255546456->zs ::radec() :rr->rr1 :th->th1 :ra->ra1 :dc->dc1 :r0->r01 :xx3->xx31 :yy3->yy31 :zz3->zz31 :xs->xs1 :ys->ys1 :zs->zs1 :e1->e11 :EndProgm rr1 = 0.550730015025 AU th1 = 1.94371466732 radians =111.366647015 deg xx31= -0.266166523942 AU yy31= 0.4618979603 AU zz31= -0.138236049035 AU ra1 = 142.094574218 degrees dc1 = 2.06588785641 degrees r01 = 1.47077284511 AU ------- :ra22() time 1=40 days :Prgm :ClIO :40->dt :dt->dt2 :-0.931108260968-xs :0.371439715781->ys :0.161052202235->zs :radec() :rr->rr2 :th->th2 :ra->ra2 :dc->dc2 :r0->r02 :xx3->xx32 :yy3->yy32 :zz3->zz32 :xs->xs2 :ys->ys2 :zs->zs2 :e1->e12 :EndProgm rr2 = 0.572141626031 AU th2 = 2.13864324044 radians = 122.535231561 degrees xx32= -0.270986191631 AU yy32= 0.439148484293 AU zz32= -0.247105509704 AU ra2 = 146.007690781 degrees dc2 = -3,39669089616 degrees r02 = 1.45240816398 AU ---- :ra33() time 3=45 days :Prgm :ClIO :45->dt :dt->dt3 :-0.962082097969-xs :0.299156786197->ys :0.129711113762->zs :radec() :rr->rr3 :th->th3 :ra->ra3 :dc->dc3 :r0->r03 :xx3->xx33 :yy3->yy33 :zz3->zz33 :xs->xs3 :ys->ys3 :zs->zs3 :e1->e13 :EndProgm rr3 = 0.591138989727 AU th3 = 132.938952042 degrees xx33= -0.265112136054 AU yy33= 0.399047678984 AU zz33= -0.3462972861445 AU ra3 = 150.362559348 degrees dc3 = -8.72114488048 degrees r03 = 1.42842736367 AU ----- Give time 1 and time 3 values, find elements Needs dt1 dt3 :xyzt() :Prgm :ClIO :SetMode("Exact/Approx","Approximate") :ra0() true elements :ra11() time 1 :ra33() time 3 :dt3-dt1->ddt :sqrt(xx31*xx31+yy31*yy31+zz31*zz31)-rr1 0.5507 :sqrt(xx33*xx33+yy33*yy33+zz33*zz33)-rr3 0.5911 :(xx31*xx33+yy31*yy33+zz31*zz33)/(rr1*rr3)->s8 0.9299 :SetMode(Angle","Radian") :0.5*cos-1(s8)->f 0.18825 :kk^2*ddt^2/(2*sqrt(rr1*rr3)*cos(f(^3->m1 0.021006 :(rr1*rr3)/(4*sqrt(rr1*rr3)*cos(f))-.5->l1 0.009311 :Define g(yb)=2*sin-1(sqrt(m1/yb^2-l1)) :zeros(yb^3-yb^2-m1*((2*g(yb)-sin(2*g(yb)))/(sin(g(yb)))^3),yb=1)->xx :product(xx)->yb :(yb*rr1*rr3*sin(2*f)/(kk*ddt))^2->pp semilatusrectum 0.5106 :Disp "f",f 0.18825 :Disp "yb" pp",{yb,pp} yb is ratio of area of segment to area of triangle :pp/rr1-1->qq1 -0.072866 :pp/rr3-1->qq3 -0.1362 :(qq1*cos(2*f)-qq3)/sin(2*f)->qq2 0.18625 :R|>Ptheta(qq1,qq2)->th1 theta 1 1.9437 radians :R|>Pr(qq1,qq2)->ee eccentricity=0.2 :pp/(1+ee)->qq perihelion = 0.4255 AU :qq/(1-ee)->aa semimajor axis AU 0.531875 :ee+rr1*cos(th1)/aa->s1 -0.17725 :rr1*sin(th1)/(aa*sqrt(1-ee*ee))->s2 :R|>Pthetaa(s1,s2)->e1 eccentric anomaly 1.748988 :(e1-ee*sin(e1))*aa^(1.5)/kk->t1 time 1 =35 :Disp "ee",ee :Disp "e1",e1 :EndPrgm ======== Find pq xyz needed to find ww,om,ii in woi() Needs xx11 yy11 zz11 xx13 yy13 zz13 th1 th3 :ab() :Prgm :ClIO :xyz() :norm([xx31,yy31,zz31})->rr1 :norm([xx33,yy33,zz33})->rr3 :Setmode("Angle","Radian") :rr1*rr3*sin(th3-th1)->s1 :(xx31*rr3*sin(th3)-xx33*rr1*sin(th1))/s1->px :(xx33*rr1*cos(th1)-xx31*rr3*cos(th3))/s1->qx :(yy31*rr3*sin(th3)-yy33*rr1*sin(th1))/s1->py :(yy33*rr1*cos(th1)-yy31*rr3*cos(th3))/s1->qy :(zz31*rr3*sin(th3)-zz33*rr1*sin(th1))/s1->pz :(zz33*rr1*cos(th1)-zz31*rr3*cos(th3))/s1->qz :EndPrgm px= 0.173630530854 qx=-0.451038790793 py=-0.039858641561 qy= 0.885007671583 pz= 0.984003926541 qz= 0.11543582823 ===== Find xyz13 Needs xs1 xs2 xs3 ys1 ys2 ys3 zs1 zs2 zs3 ro1 ro2 ro3 :xyz() :Prgm :ClIO :Setmode("Angle","Degree") :abc() :r01*aa1-xs1->xx1 :r01*bb1-ys1->yy1 :r01*cc1-ys1->zz1 :r02*aa2-xs2->xx2 :r02*bb2-ys2->yy2 :r02*cc2-ys2->zz2 :r03*aa3-xs3->xx3 :r03*bb3-ys3->yy3 :r03*cc3-ys3->zz3 :EndPrgm ===== Find abc direction cosines Needs ra1 dc1 ra2 dc2 ra3 dc3 :abc() :Prgm :ClIO :cos(ra1)*cos(dc1)->aa1 :sin(ra1)*cos(dc1)->bb1 :sin(dc1)->cc1 :cos(ra2)*cos(dc2)->aa2 :sin(ra2)*cos(dc2)->bb2 :sin(dc2)->cc2 :cos(ra3)*cos(dc3)->aa3 :sin(ra3)*cos(dc3)->bb3 :sin(dc3)->cc3 :EndPrgm ====== Find ww,om, ii from pqxyz from ab() ie :woi() :Prgm :ab() :ClIO :Setmode("Angle","Degree") :pz*cos(ie)-py*sin(ie)->s1 :qz*cos(ie)-qy*sin(ie)->s2 :R|>Pr(s1,s2)->s3 :R|>Ptheta(s2,s1)->ww :(py*cos(ww)-qy*sin(ww))/(cos(ie))->s4 :(px*cos(ww)-qx*sin(ww))/(cos(ie))->s5 :R|>Ptheta(s5,s4)->om :-(px*sin(ww)+qx*cos(ww))/sin(om)->s6 :R|>Ptheta(s6,s3)->ii :Disp "ww",ww :Disp "om",om :Disp "ii",ii :EndPrgm ------------ s1=.91865005135 s2=-.246151539386 s3=.9510565516295 s4=-.920504853452 s5=.390731128489 Elements ww=105. om=-67. ii=72. ===== Given true elements, find velocity and position. Given velocity and position., find true elements. sqrt is square root symbol cos-1() is arccos |> is right arrowhead theta is theta symbol Program elv(0 needs vel(), ra0(), ra22(), woi() and radec() :elv() :Prgm :ClrIO :SetMode("Exact/Approx","Approximate") :vel() :[[xx3,yy3,zz3]]->s1 :[[dx3,dy3,dz3]]->s2 :Disp "s1,s2",s1,s2 :Pause :ClrIO :norm([xx3,yy3,zz3])->rr2 :norm([dx3,dy3,dz3])->vl :Pause :ClrIO :-kk^2/(vl^2-2*kk^2/rr2)->aa :dotP(s1,s2)->s3 :Disp "aa,s3",aa,s3 :Pause :ClrIO :s3/(rr2*abs(vl))->s3a :cos-1(s3a)->s4 :Disp "s3a,s4",s3a,s4 :Pause :ClrIO :vl*sin(s4)/rr2->dth :vl*cos(s4)->drt :Disp "dth,drt",dth,drt :Pause :ClrIO :(dth*rr2^2/kk)^2->pp :pp/rr2-1->s5 :Disp "pp,s5",pp,s5 :Pause :ClrIO :drt*sqrt(pp)/kk->s6 :R|>Pr(s5,s6)->ee :Disp "s6,ee",s6,ee :Pause :ClrIO :R|>Ptheta(s5,s6)->th :rr2^2*cos(th)/aa+ee->s7 :Disp "th,s7",th,s7 :Pause :ClrIO :rr2*sin(th)/(aa*sqrt(1+ee^2))->s8 :SetMode("Angle","Radian") :R|>Ptheta(s7,s8)->e1 :Disp "s8,e1",s8,e1 :Pause :ClrIO :(e1-ee*sin(e1))*aa^(1.5)/kk->tt :SetMode("Angle","Degree") :rr2*cos(th)->s1 :Disp "tt,s1",tt,s1 :Pause :ClrIO :rr2*sin(th)->s2 :drt*cos(th)-rr2*sin(th)*dth->s3 :Disp "s2,s3",s2,s3 :Pause :ClrIO :drt*sin(th)+rr2*cos(th)*dth->s4 :Disp "s4",s4 :Pause :ClrIO :[[s1,s2][s3,s4]]->s6 :s6^(-1)->s7 :Disp "s6,s7",s6,s7 :Pause :ClrIO :[[xx3][dx3]]->s8 :s7*s8->s9 :Disp "s8,s9",s8,s9 :Pause :ClrIO :subMat(s9,1,1,1,1)->px :subMat(s9,2,1,2,1)->qx :det(px)->px :det(qx)->qx :Disp "px,qx",px,qx :Pause :ClrIO :[[yy3][dy3]]->s8 :s7*s8->s9 :Disp "s9",s9 :Pause :ClrIO :subMat(s9,1,1,1,1)->py :subMat(s9,2,1,2,1)->qy :Disp "py,qy",py,qy :Pause :ClrIO :[[zz3][dz3]]->s8 :s7*s8->s9 :Disp "s9",s9 :Pause :ClrIO :subMat(s9,1,1,1,1)->pz :subMat(s9,2,1,2,1)->qz :Disp "pz,qz",pz,qz :Pause :ClrIO :woi() :Disp "tt,ee",tt,ee :Disp "qq",qq :EndPrgm ---- variables generated by elv() aa= 0.531875 semimajor axis dc=dc2= -3.39669019616 dddy= {1.95900897925} eccentric anomaly E drt= 0.004059102127 = AU/day dt=dt2= 40 days after perihelion dth=0.037550416193 = d(th)/dt dx1= 0.013108254602 dx3= 0.000144151811 dy= -0.0174991124294 dy1= -0.005407523147 dy2= -0.014179099603 dy3= -0.006388371699 dz2= -0.016642644967 dz3= -0.020909643806 e1= e12 = 1.95900897925 eccentric anomaly E ee= 0.2 eccentricity ie=23.441028 deg = tilt of earth ii= 72 deg = inclination of orbit kk= 0.01720209895 Gauss gravity constant mm= 1.77389155705 Mean anomaly radians om= -67 Longitude of ascending node degrees pp= 0.5106 AU semilatus rectum px= 0.173630530854 py=-0.039858641561 pz=0.984003926541 qq= 0.4255 = perihelion AU qx= -0.451038790793 qy= 0.885007671583 qz= 0.11543582823 r0= r02= 1.45240816398 ra= ra2= 146.007690781 Right Ascension degrees ram= 0.503897334963 rom= 0.40769864625 rr= rr2 =0.572141626031 radius s1= 0.91865005135 s2= -0.246151539386 s3= 0.951056516295 s3a= 0.185650217548 s4= -0.920504853452 s5= 0.390731128489 s6= 0.309016994375 ss1= -1.2020944526 ss2= 0.810588200074 ss3= -0.086653307469 th=th2= 122.535731561 degrees true anomaly tha= -52.8071544636 tt=40 days past perihelion u=226.535231561 vl= 0.02186424654 ww= 105 degrees = argument of perigee x= -0.307708130153 x1= -0.386273822554 xs= xs2= -0.931108260968 x earth to sun AU xx1= -0.386273822554 xx3= xx32= -0.270986191631 y= -0.422064656477 ys= ys2= 0.371439715781 y earth to sun AU yy0= -0.422064656477 yy1= -0.130425151576 yy2= 0.304605761763 yy3= yy32= 0.439148484293 zs=zs2= 0.161052202235 z earth to sun AU zz= zz2= -0.40140734184 zz3= zz32= -0.247105509704 ======== Find xyz velocity time 40 days :vel() :Prgm :ClrIO :ra0() :ra22() :qq*(1+ee)->pp :kk*sqrt(pp)/rr2^2->dth :kk*ee*sin(th2)/(sqrt(pp))->drt :th2+ww->u :rr2*cos(u)->x1 :rr2*sin(u)->y :drt*cos(u)-rr2*sin(u)*dth->dx1 :drt*sin(u)*rr2*cos(u)*dth->dy :y*cos(ii)->yy1 :y*sin(ii)->zz2 :dy*sin(ii)->dz2 :dy*cos(ii)->dy1 :dx1*cos(om)-dy1*sin(om)->dx3 :dy1*cos(om)+dx1*sin(om)->dy2 :dy2*cos(ie)-dz2*sin(ie)->dy3 :dz2*cos(ie)+dy2*sin(ie)->dz3 :EndPrgm ======= Run programs in this order ra0() True elemwnts ra11() ra, dc r theta xyz3 ra22() ra33() ab() Find pqxyz ww,om,ii xyzt() find e1, t1, aa, qq ====