 | |
| Red circle - Mars impact locus. Green star - old nominal approach. Black ellipses - old 1,2,3 sigma. Blue dots - new Monte Carlo samples. Orange rings - new 1,2,3 sigma ellipses according to Monte Carlo |
Exact steps to duplicate what I am doing, for the benefit of posterity.
- Get Find_orb from Project Pluto. I got the Win32 version in the Download section of that page.
- Get observations from the Minor Planet Center. It seems like they only update once per week. There is a download observations link at the top of the observations table at the bottom of the page. Get that file.
- Append the observations file to the example file in the Find_orb package
- Run find_o32.exe
- Click open and open the file example . The list of objects should include C/2013 A1.
- Double click on C/2013 A1, and the machine will calculate a preliminary orbit centered on an epoch in 2012 (for this object, perhaps a different epoch for a different object)
- Change the epoch to 2014 Oct 19.8 (or whenever the close approach is)
- Turn on all the perturber check boxes
- Click the Monte Carlo button, then wait for it to accumulate about 5000 samples. It calculates a few orbits per second, so this may take 10-20 minutes.
- Results are saved in a file state.txt . Run that through a pair of filters to get the position and velocity vectors into separate files pos.txt and vel.txt .
grep "#.* AU" state.txt > pos.txt
grep "#.* mAU" state.txt > vel.txt
- Run pos.txt and vel.txt through the following IDL script. Getting IDL is beyond the scope of this document, but if you don't want to shell out real cash for IDL, consider GDL, gnu data language.
- Comparison is with JPL small body database.
IDL script used:
function asinh,z
return,alog(z+sqrt(1d +z^2))
end
;Calculate impact parameter B from gravity parameter of planet and
;entry interface conditions
;input
; mu - gravity parameter of target planet. Implies units of distance and time
; r_ei - Entry interface distance from center of planet in distance units implied by mu
; fpa_ei - inertial flight path angle at entry interface in radians, negative is downward (usual case)
; v_ei= - inertial speed at entry interface in distance and time units implied by mu
; v_inf= - hyperbolic excess speed, in distance and time units implied by mu
; One or the other of v_ei= and v_inf= must be passed as input, the other is calculated and returned as output
; if both are passed in, v_inf is recalculated from v_ei
;output
; v_ei= or v_inf= - one must be passed in, the other is calculated and passed out
; ta_inf - true anomaly at infinity in radians, always negative
; ta_ei - true anomaly at entry interface altitude in radians, always negative
; dt_ei_rp - time from entry interface to theoretical periapse, ignoring aerobraking and lithobraking
;return
; impact parameter B - distance from center of planet to hyperbolic asymtote in
; distance units implied by mu
function b_from_ei,mu,r_ei,fpa_ei,v_ei=v_ei,v_inf=v_inf,ta_inf=ta_inf,ta_ei=ta_ei,dt_ei_rp=dt_ei_rp,rp=rp
tau=!dpi*2 ;Tau manifesto
if n_elements(v_ei) ne 0 then begin
xi=v_ei^2/2-mu/r_ei
a=-mu/(2*xi)
v_inf=sqrt(-mu/a)
end else begin
; specific energy at infinity, remembering that potential energy is
; zero at infinity and is negative at finite distance
xi=v_inf^2/2
a=-mu/(2*xi)
v_ei=sqrt(2*mu/r_ei-mu/a)
end
h=r_ei*v_ei*sin(tau/4d -fpa_ei) ;magnitude of cross product a x b is a*b*sin(theta)
p=h^2/mu
e=sqrt(1d +(2d*xi*h^2)/(mu^2))
rp=a*(1d -e)
B=rp*sqrt((2d*mu)/(rp*v_inf^2)+1d)
ta_ei=-acos(p/(r_ei*e)-1d/e)
ta_inf=-acos(-1/e)
shEE=sin(ta_ei)*sqrt(e^2-1)/(1d +e*cos(ta_ei))
MM=e*shEE-asinh(shEE)
n=sqrt(-mu/(a^3))
dt_ei_rp=MM/n
return,b
end
function b_plane,rv,vv,l_du,mu
tau=!dpi*2d ; tau manifesto
r=vlength(rv)
v=vlength(vv)
fpa=tau/4d -acos(dotp(rv,vv)/(r*v))
b=b_from_ei(mu,r,fpa,v_ei=v,v_inf=v_inf)
if n_elements(l_du) eq 0 then begin
b=su_to_cu(b,l_du,mu,1,0,/inv)
end
return,b
end
;Component magnitude of vector u in the direction of vector b
function comp,u,b
return,dotp(u,b)/vlength(b)
end
;A Grid is a multidimensional array of vectors. It is inconvenient to talk about
;an array of things which themselves are arrays, so it'ss a grid.
;
;A grid of vectors can have any number of dimensions, but the last dimension
;is the vector component index. So, a grid would typically look like v[i,j,3].
;Even though this is a 3D array, it is only a 2D grid of vectors.
;
;A grid can have any number of dimensions, but if more than three are needed, this function
;and resolve_grid will need to be extended.
;
;This function takes three arrays of scalars of the same size and shape, and builds a grid
;of vectors with the same size and shape. Each input array becomes one vector component.
function compose_grid,x,y,z,w,v,u,t,s
ndims=size(x,/n_dimensions)
ss=0
if n_elements(x) gt 0 then ss+=1
if n_elements(y) gt 0 then ss+=1
if n_elements(z) gt 0 then ss+=1
if n_elements(w) gt 0 then ss+=1
if n_elements(v) gt 0 then ss+=1
if n_elements(u) gt 0 then ss+=1
if n_elements(t) gt 0 then ss+=1
if(ndims gt 0) then begin
result=make_array([size(x,/dimensions),ss],type=size(x,/type));
end else begin
result=make_array(ss,type=size(x,/type));
end
case ndims of
0: begin
if ss ge 1 then result[0]=x;
if ss ge 2 then result[1]=y;
if ss ge 3 then result[2]=z;
if ss ge 4 then result[3]=w;
if ss ge 5 then result[4]=v;
if ss ge 6 then result[5]=u;
if ss ge 7 then result[6]=t;
if ss ge 8 then result[7]=s;
end
1: begin
if ss ge 1 then result[*,0]=x;
if ss ge 2 then result[*,1]=y;
if ss ge 3 then result[*,2]=z;
if ss ge 4 then result[*,3]=w;
if ss ge 5 then result[*,4]=v;
if ss ge 6 then result[*,5]=u;
if ss ge 7 then result[*,6]=t;
if ss ge 8 then result[*,7]=s;
end
2: begin
if ss ge 1 then result[*,*,0]=x;
if ss ge 2 then result[*,*,1]=y;
if ss ge 3 then result[*,*,2]=z;
if ss ge 4 then result[*,*,3]=w;
if ss ge 5 then result[*,*,4]=v;
if ss ge 6 then result[*,*,5]=u;
if ss ge 7 then result[*,*,6]=t;
if ss ge 8 then result[*,*,7]=s;
end
3: begin
if ss ge 1 then result[*,*,*,0]=x;
if ss ge 2 then result[*,*,*,1]=y;
if ss ge 3 then result[*,*,*,2]=z;
if ss ge 4 then result[*,*,*,3]=w;
if ss ge 5 then result[*,*,*,4]=v;
if ss ge 6 then result[*,*,*,5]=u;
if ss ge 7 then result[*,*,*,6]=t;
if ss ge 8 then result[*,*,*,7]=s;
end
4: begin
if ss ge 1 then result[*,*,*,*,0]=x;
if ss ge 2 then result[*,*,*,*,1]=y;
if ss ge 3 then result[*,*,*,*,2]=z;
if ss ge 4 then result[*,*,*,*,3]=w;
if ss ge 5 then result[*,*,*,*,4]=v;
if ss ge 6 then result[*,*,*,*,5]=u;
if ss ge 7 then result[*,*,*,*,6]=t;
if ss ge 8 then result[*,*,*,*,7]=s;
end
5: begin
if ss ge 1 then result[*,*,*,*,*,0]=x;
if ss ge 2 then result[*,*,*,*,*,1]=y;
if ss ge 3 then result[*,*,*,*,*,2]=z;
if ss ge 4 then result[*,*,*,*,*,3]=w;
if ss ge 5 then result[*,*,*,*,*,4]=v;
if ss ge 6 then result[*,*,*,*,*,5]=u;
if ss ge 7 then result[*,*,*,*,*,6]=t;
if ss ge 8 then result[*,*,*,*,*,7]=s;
end
6: begin
if ss ge 1 then result[*,*,*,*,*,*,0]=x;
if ss ge 2 then result[*,*,*,*,*,*,1]=y;
if ss ge 3 then result[*,*,*,*,*,*,2]=z;
if ss ge 4 then result[*,*,*,*,*,*,3]=w;
if ss ge 5 then result[*,*,*,*,*,*,4]=v;
if ss ge 6 then result[*,*,*,*,*,*,5]=u;
if ss ge 7 then result[*,*,*,*,*,*,6]=t;
if ss ge 8 then result[*,*,*,*,*,*,7]=s;
end
7: begin
if ss ge 1 then result[*,*,*,*,*,*,*,0]=x;
if ss ge 2 then result[*,*,*,*,*,*,*,1]=y;
if ss ge 3 then result[*,*,*,*,*,*,*,2]=z;
if ss ge 4 then result[*,*,*,*,*,*,*,3]=w;
if ss ge 5 then result[*,*,*,*,*,*,*,4]=v;
if ss ge 6 then result[*,*,*,*,*,*,*,5]=u;
if ss ge 7 then result[*,*,*,*,*,*,*,6]=t;
if ss ge 8 then result[*,*,*,*,*,*,*,7]=s;
end
else: begin
print,"Unsupported grid dimension"
return,-1
end
endcase
return,result
end
Function crossp_grid,a,b
;
;+
; NAME:
; CROSSP
;
; PURPOSE:
; Evaluates the cross product of v1 with v2
;
; CATEGORY:
; Vector mathematics.
;
; CALLING SEQUENCE:
; Result = CROSSP(v1, v2)
;
; INPUTS:
; v1, v2: Three-element vectors.
;
; OUTPUTS:
; Returns a 3-element, floating-point vector.
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
; None.
;
; RESTRICTIONS:
; Vectors must have 3 elements.
;
; PROCEDURE:
; v1 X v2 = | i j k | = (b1c2 - b2c1)i + (c1a2-c2a1)j + (a1b2-a2b1)k
; | a1 b1 c1 |
; | a2 b2 c2 |
;
; MODIFICATION HISTORY:
; Written, DMS, Aug, 1983;
; Modified by Chris Jeppesen, Jul, 2005 to multiply a whole grid
;-
resolve_grid,a,x=ax,y=ay,z=az
resolve_grid,b,x=bx,y=by,z=bz
x=ay*bz-az*by
y=az*bx-ax*bz
z=ax*by-ay*bx
return,compose_grid(x,y,z)
end
function dotp,a,b,dim
resolve_grid,a,x=ax,y=ay,z=az,w=aw,v=av,u=au,t=at,s=as
resolve_grid,b,x=bx,y=by,z=bz,w=bw,v=bv,u=bu,t=bt,s=bs
acc=ax*bx
if n_elements(ay) gt 0 and n_elements(by) gt 0 then acc+=ay*by
if n_elements(az) gt 0 and n_elements(bz) gt 0 then acc+=az*bz
if n_elements(aw) gt 0 and n_elements(bw) gt 0 then acc+=aw*bw
if n_elements(av) gt 0 and n_elements(bv) gt 0 then acc+=av*bv
if n_elements(au) gt 0 and n_elements(bu) gt 0 then acc+=au*bu
if n_elements(at) gt 0 and n_elements(bt) gt 0 then acc+=at*bt
if n_elements(as) gt 0 and n_elements(bs) gt 0 then acc+=as*bs
return,acc
end
;Given state vector, calculate orbital elements
;input
; r - position vector, may be a grid of vectors, in distance units implied by mu, any inertial frame is fine, but center of attraction is assumed to be
; at the origin of the frame
; v - inertial velocity vector, must follow r as either not a grid or a grid of the same size, distance and time units implied by mu
; must be in same frame as r
; l_du - length of a distance unit, used for conversion to canonical units internally
; mu - gravity parameter, implies distance and time units
;output
; ev= eccentricity vector. Follows grid-ness of r. Length is equal to eccentricity, points from center of attraction to periapse
;return
; a structure, fields of which are either scalar or grids depending on grid-ness of r
; p: semi-parameter, distance from focus to orbit at TA=+-90deg, in original distance units, always positive for any eccentricity
; a: semimajor axis, in original distance units
; e: eccentricity
; i: inclination, radians
; an: longitude of ascending node, angle between x axis and line of intersection between orbit plane and xy plane, radians
; ap: argument of periapse, angle between xy plane and periapse along orbit plane, radians
; ta: true anomaly, angle between periapse and object, radians
; tp: time to next periapse original time units. Negative if only one periapse and in the past
; rp: radius of periapse in original distance units
function elorb,r_,v_,l_DU,mu,ev=ev
tau=2d*!dpi ;Tau manifesto
rv=su_to_cu(r_,l_DU,mu,1,0)
vv=su_to_cu(v_,l_DU,mu,1,-1)
r=vlength(rv)
v=vlength(vv)
hv=crossp_grid(rv,vv)
h=vlength(hv)
nv=crossp_grid([0,0,1],hv)
n=vlength(nv)
ev=smult_grid((v^2-1d/r),rv)-smult_grid(dotp(rv,vv),vv)
e=vlength(ev)
xi=v^2/2d -1/r
if n_elements(xi) eq 1 then begin
a=0d
p=0d
end else begin
a=dblarr(size(xi,/dim))
p=dblarr(size(xi,/dim))
end
w=where(e ne 1.0,count)
if count gt 0 then begin
a[w]=-1/(2*xi[w])
p[w]=a[w]*(1-e[w]^2)
end
w=where(e eq 1.0,count)
if count gt 0 then begin
;parabolic case
p[w]=h[w]^2
a[w]=!values.d_infinity
end
resolve_grid,hv,x=hx,y=hy,z=hz
resolve_grid,nv,x=nx,y=ny,z=nz
i=acos(hz/h)
an=acos(nx/n)
w=where(ny lt 0,count)
if count gt 0 then an[w]=tau-an[w]
ap=vangle(ev,nv)
resolve_grid,ev,x=ex,y=ey,z=ez
w=where(ez lt 0,count)
if count gt 0 then ap[w]=tau-ap[w]
ta=vangle(ev,rv)
w=where(dotp(rv,vv) lt 0,count)
if count gt 0 then ta[w]=tau-ta[w]
if n_elements(a) eq 1 then begin
rp=0d
n=0d
MM=0d
end else begin
rp=dblarr(size(a,/dim))
n=dblarr(size(a,/dim))
MM=dblarr(size(a,/dim))
end
w=where(a gt 0,count)
if count gt 0 then begin
EE=2*atan(sqrt((1d -e[w])/(1d +e[w])*tan(ta[w]/2d)))
MM[w]=EE-e[w]*sin(EE)
n[w]=sqrt(1d/(a[w]^3))
rp[w]=a[w]*(1d -e)
end
w=where(a lt 0,count)
if count gt 0 then begin
EE=asinh(sin(ta[w])*sqrt(e[w]^2-1)/(1d +e[w]*cos(ta[w])))
MM[w]=e[w]*sinh(EE)-EE
n[w]=sqrt(-1d/(a[w]^3))
rp[w]=a[w]*(1d -e[w])
end
w=where(~finite(a),count)
if count gt 0 then begin
EE=tan/2
MM[w]=EE^3/3d +EE
n[w]=sqrt(2d/(p[w]^3))
rp[w]=p[w]/2
end
tp=MM/n
return,{p:su_to_cu(p,l_du,mu,1,0,/inv), $
a:su_to_cu(a,l_du,mu,1,0,/inv), $
e:e, $
i:i, $
an:an, $
ap:ap, $
ta:ta, $
tp:su_to_cu(tp,l_du,mu,0,1,/inv), $
rp:su_to_cu(rp,l_du,mu,1,0,/inv) $
}
end
function normalize_grid,v
resolve_grid,v,x=x,y=y,z=z
len=vlength(v);
result=compose_grid(x/len,y/len,z/len)
return,result;
end
;Orthogonal projection of vector u onto vector b, or vector component of
;u in the direction of b
function proj,u,b
return,smult_grid(comp(u,b)/vlength(b),b)
end
;Orthogonal projection of vector u not onto vector b, or vector component of
;u perpendicular to the direction of b
function proj_p,u,b
return,u-proj(u,b)
end
function smult_grid,s,v
if n_elements(s) eq 1 then return,s*v
return,rebin(s,size(v,/dim))*v
end
pro resolve_grid,vv,x=x,y=y,z=z,w=w,v=v,u=u,t=t,s=s,n_dimension_vec=n_dimension_vec
ss=size(vv,/dim)
ndims=n_elements(ss)-1
n_dimension_vec=ss[ndims]
;Column vector exception
if ndims eq 1 and ss[0] eq 1 then ndims--
case ndims of
-1: begin
message,"Must pass a vector or grid of vectors"
return
end
0: begin
if n_dimension_vec ge 1 then x=vv[0];
if n_dimension_vec ge 2 then y=vv[1];
if n_dimension_vec ge 3 then z=vv[2];
if n_dimension_vec ge 4 then w=vv[3];
if n_dimension_vec ge 5 then v=vv[4];
if n_dimension_vec ge 6 then u=vv[5];
if n_dimension_vec ge 7 then t=vv[6];
if n_dimension_vec ge 8 then s=vv[7];
end
1: begin
if n_dimension_vec ge 1 then x=vv[*,0];
if n_dimension_vec ge 2 then y=vv[*,1];
if n_dimension_vec ge 3 then z=vv[*,2];
if n_dimension_vec ge 4 then w=vv[*,3];
if n_dimension_vec ge 5 then v=vv[*,4];
if n_dimension_vec ge 6 then u=vv[*,5];
if n_dimension_vec ge 7 then t=vv[*,6];
if n_dimension_vec ge 8 then s=vv[*,7];
end
2: begin
if n_dimension_vec ge 1 then x=vv[*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,7];
end
3: begin
if n_dimension_vec ge 1 then x=vv[*,*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,*,7];
end
4: begin
if n_dimension_vec ge 1 then x=vv[*,*,*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,*,*,7];
end
5: begin
if n_dimension_vec ge 1 then x=vv[*,*,*,*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,*,*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,*,*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,*,*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,*,*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,*,*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,*,*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,*,*,*,7];
end
6: begin
if n_dimension_vec ge 1 then x=vv[*,*,*,*,*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,*,*,*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,*,*,*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,*,*,*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,*,*,*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,*,*,*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,*,*,*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,*,*,*,*,7];
end
7: begin
if n_dimension_vec ge 1 then x=vv[*,*,*,*,*,*,*,0];
if n_dimension_vec ge 2 then y=vv[*,*,*,*,*,*,*,1];
if n_dimension_vec ge 3 then z=vv[*,*,*,*,*,*,*,2];
if n_dimension_vec ge 4 then w=vv[*,*,*,*,*,*,*,3];
if n_dimension_vec ge 5 then v=vv[*,*,*,*,*,*,*,4];
if n_dimension_vec ge 6 then u=vv[*,*,*,*,*,*,*,5];
if n_dimension_vec ge 7 then t=vv[*,*,*,*,*,*,*,6];
if n_dimension_vec ge 8 then s=vv[*,*,*,*,*,*,*,7];
end
else: begin
message,"Unsupported grid dimension"
end
endcase
end
;Convert standard units to canonical units
;input
; x - measurement to convert, may be scalar or any size or shape of vector
; a - length of canonical distance unit in standard units. Standard length unit is implied by this.
; mu - gravitational constant in standard distance and time units. Standard length unit same as
; above, standard time unit implied by this.
; LL - power of length dimension used in x
; TT - power of time dimension used in x
; /inverse - convert x from canonical units back to standard units
;return
; a scalar or array the same size as x, in canonical units (or standard if /inv was set)
;Example
; An object orbiting Earth has a position of <1131340,-2282343,6672423> m
; and a speed of <-5643.05,4303.33,2428.79> m/s. Convert this to canonical units
; Earth radius used as distance unit length: 6378137m
; Earth gravitational constant: 398600.4415d9 m,s
; print,su_to_cu([1131340d,-2282343d,6672423d],6378137d,398600.4415d9,1,0)
; 0.17737781 -0.35783850 1.0461398
; print,su_to_cu([-5643.05d,4303.33d,2428.79d],6378137d,398600.4415d9,1,-1)
; -0.71382529 0.54435559 0.30723310
; We are going to solve the Kepler problem over a time of 40min=2400s. How many canonical time units?
; print,su_to_cu(2400d,6378137d,398600.4415d9,0,1)
; 2.9746739
; The answer is r_t=<-0.6616125, 0.6840739,-0.6206809> and
; v_t=< 0.4667380,-0.2424455,-0.7732126>. What is this in SI units?
; print,su_to_cu([-0.6616125d, 0.6840739d,-0.6206809d],6378137d,398600.4415d9,1,0,/inv)
; -4219855.2 4363117.1 -3958787.8
; print,su_to_cu([ 0.4667380d,-0.2424455d,-0.7732126d],6378137d,398600.4415d9,1,-1,/inv)
; 3689.7346 -1916.6203 -6112.5284
function su_to_cu,x,a,mu,LL,TT,inverse=inverse
;To convert: Standard Unit to Canonical Unit Multiply by
;Distance m DU 1/a
;Time s TU sqrt(mu/a^3)
;For derived units, raise the base unit for each dimension to the power of the dimension needed, then multiply. For example
;Speed m/s DU/TU (LL=1,TT=-1) 1/a*sqrt(a^3/mu)
;For inverse conversion, divide instead of multiply
DU=(1d/a[0])^LL
TU=sqrt(mu[0]/a[0]^3)^TT
DUTU=DU*TU
if keyword_set(inverse) then DUTU=1d/DUTU
return,x*DUTU
end
;Function returns the angle between two vectors in radians.
function vangle,a,b,dim
return,acos(dotp(a,b,dim)/(vlength(a,dim)*vlength(b,dim)))
end
function vlength,v,dim
return,sqrt(dotp(v,v,dim))
end
pro read_monte_carlo
template={ $
VERSION: 1.0000000e+000, $
DATASTART: 0L, $
DELIMITER: 32B, $
MISSINGVALUE: !values.f_nan, $
COMMENTSYMBOL: '', $
FIELDCOUNT: 5L, $
FIELDTYPES: [0L,5L,5L,5L,0L], $
FIELDNAMES: ['FIELD1','x','y','z','FIELD5'], $
FIELDLOCATIONS: [0L,5L,21L,38L,54L], $
FIELDGROUPS: [0L,1L,2L,3L,4L] $
}
pos=read_ascii('pos.txt',template=template)
vel=read_ascii('vel.txt',template=template)
AU=149597870.691d
D=86400d
state=[[pos.x*AU], $
[pos.y*AU], $
[pos.z*AU], $
[vel.x/1000d*AU/D], $
[vel.y/1000d*AU/D], $
[vel.z/1000d*AU/D]]
help,state
cspice_gdpool,'BODY10_GM',0,1,sun_mu,found
sun_mu=sun_mu[0]
cspice_gdpool,'BODY499_GM',0,1,mars_mu,found
mars_mu=mars_mu[0]
cspice_gdpool,'BODY499_RADII',0,1,mars_r,found
mars_r=mars_r[0]
el_sun=elorb(state[*,0:2],state[*,3:5],AU,sun_mu,ev=ev)
start_cal='2456950.3JDTDT'
cspice_str2et,start_cal,et
cspice_spkezr,'499',et,'ECLIPJ2000','NONE','10',mars_state,ltime
state=state-rebin(transpose(mars_state),n_elements(pos.x),6)
el_mars=elorb(state[*,0:2],state[*,3:5],mars_r,mars_mu,ev=ev)
peri=smult_grid(1d/vlength(ev)*el_mars.rp,ev)
sh=normalize_grid(state[*,3:5])
th=normalize_grid(crossp_grid(sh,[0,0,-1]))
rh=normalize_grid(crossp_grid(th,sh))
b=b_plane(state[*,0:2],state[*,3:5],mars_r,mars_mu)
bv=smult_grid(b,normalize_grid(proj_p(peri,sh)))
bdotr=dotp(bv,rh)
bdott=dotp(bv,th)
; bdotr=bdotr[0:100]
; bdott=bdott[0:100]
thetab=atan(bdotr,bdott)
s=3d5
plot,bdott,bdotr,psym=3,color=64,xrange=[-s,s],yrange=[s,-s],/iso,/xs,/ys,xtitle="B.t km",ytitle="B.r km",charsize=2,title="Monte Carlo uncertainty"
;oplot,bdott,bdotr,psym=3,color=64;,xrange=[-s,s],yrange=[s,-s];,/iso,xtitle="B.t km",ytitle="B.r km",psym=3,charsize=2,title="Monte Carlo uncertainty
q=dindgen(!dpi*200)/100d
; oplot,cos(q)*mars_r,sin(q)*mars_r
c=cos(q)
s=sin(q)
f=[[bdott],[bdotr]]
xbar=total(f,1)/double(n_elements(bdotr))
help,xbar
print,xbar
oplot,[1d,1]*xbar[0],[1d,1]*xbar[1],psym=2
fac=F - rebin(transpose(xbar),n_elements(bdotr),n_elements(xbar))
P=(Fac ## transpose(Fac))/(n_elements(bdotr)-1)
L=P
la_choldc,L
n=2
for i=0,n-2 do L[i+1:*,i]=0
print,L
Y=L ## [[c],[s]]
for i=1,3 do oplot,i*Y[*,0]+xbar[0],i*Y[*,1]+xbar[1],color=220
end
No comments:
Post a Comment