The IDL code at the end of this message when saved as a file cosmic.pro
will plot the trajectory of a cosmic ray in the earth's magnetic field in
3 dimensions. It plots the trajectory as it is happening(not at the end),
which is the way one should do during a simulation.
I would like to draw a shaded sphere (even better a sphere with earth map
on its surface) of radius rearth so that the subsequent cosmic ray trajectory
in the earth's magnetic field appears properly in relation to that sphere.
We are in the process of switching from MATLAB to IDL and I am
converting a large number of teaching programs.
I have a routine that draws the sphere (from the POLYSHADE help) when I
use it as a separate program.
Sphere.pro code:
^^^^^^^^^^^^^^^^
SPHERE = FLTARR(20, 20, 20)
FOR X=0,19 DO FOR Y=0,19 DO FOR Z=0,19 DO $
SPHERE(X, Y, Z) = SQRT((X-10)^2 + (Y-10)^2 + (Z-10)^2)
SHADE_VOLUME, SPHERE, 8, V, P
SCALE3, XRANGE=[0,20], YRANGE=[0,20], ZRANGE=[0,20]
image = POLYSHADE(V, P, /T3D) ;Render the image.
TV, image
But when I try to insert it into the cosmic.pro code
I get all kinds of conflicts between the shade_volume, scale3 and
polyshade routines and the routines I am using in cosmic.pro.
I am not experienced enough with the IDL 3-D stuff yet to figure out how
to do it.
This easy to do in MATLAB and I am sure it is just as easy to do in IDL.
Section of the old MATLAB program:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
load topo
[X,Y,Z]=sphere(24);
X=rearth*X;
Y=rearth*Y;
Z=rearth*Z;
q=surface(X,Y,Z,'FaceColor','texture','CData',topo);
rotate(q,[0 90],-45,[0 0 0]);
view(-135,25);
colormap(topomap1);
Any ideas would be helpful.
Thanks,
John Boccio
Department of Physics
Swarthmore College
boccio@swarthmore.edu
Cosmic ray trajectory code:
^^^^^^^^^^^^^^^^^^^^^^^^^^^
function accelx,alpha,r,x,y,z,vx,vy,vz
return,(alpha/r^5)*((2.0*z*z-x*x-y*y)*vy-3.0*y*z*vz)
end
function accely,alpha,r,x,y,z,vx,vy,vz
return,(alpha/r^5)*(3.0*x*z*vz-(2.0*z*z-x*x-y*y)*vx)
end
function accelz,alpha,r,x,y,z,vx,vy,vz
return,3.0*(alpha/r^5)*z*(y*vx-x*vy)
end
pro cosmic
rearth=6.4E6
x=3.0*rearth
vx=-0.9E4
y=3.0*rearth
vy=-0.9E4
z=0.0
vz=-2.9E4
h=2.0
window,0,xsize=500,ysize=500, title='Cosmic Rays'
red=[0,1,1,0,0,1]
green=[0,1,0,1,0,1]
blue=[0,1,0,0,1,0]
tvlct, 255*red,255*green,255*blue
surface,fltarr(2,2),/nodata,/save,xtitle='X', $
ytitle='Y',ztitle='Z', $
xrange=[-3.0*rearth,3.0*rearth],yrange=[-3.0*rearth,3.0*rear th], $
zrange=[-3.0*rearth,3.0*rearth],az=60.0,ax=30.0
x1=[x,x]
y1=[y,y]
z1=[z,z]
plots,x1,y1,z1,psym=3,color=2,/t3d
alpha=1.0E20
r=sqrt(x*x+y*y+z*z)
count=0L
while ((r GT rearth) AND (count LT 200000) AND (r LT 6*rearth)) do begin
count=count+1
fx1=vx
gx1=accelx(alpha,r,x,y,z,vx,vy,vz)
fy1=vy
gy1=accely(alpha,r,x,y,z,vx,vy,vz)
fz1=vz
gz1=accelz(alpha,r,x,y,z,vx,vy,vz)
fx2=vx+h*gx1/2.0
gx2=accelx(alpha,r,x+h*fx1/2.0,y+h*fy1/2.0,z+h*fz1/2.0, $
vx+h*gx1/2.0,vy+h*gy1/2.0,vz+h*gz1/2.0)
fy2=vy+h*gy1/2.0
gy2=accely(alpha,r,x+h*fx1/2.0,y+h*fy1/2.0,z+h*fz1/2.0, $
vx+h*gx1/2.0,vy+h*gy1/2.0,vz+h*gz1/2.0)
fz2=vz+h*gz1/2.0
gz2=accelz(alpha,r,x+h*fx1/2.0,y+h*fy1/2.0,z+h*fz1/2.0, $
vx+h*gx1/2.0,vy+h*gy1/2.0,vz+h*gz1/2.0)
fx3=vx+h*gx2/2.0
gx3=accelx(alpha,r,x+h*fx2/2.0,y+h*fy2/2.0,z+h*fz2/2.0, $
vx+h*gx2/2.0,vy+h*gy2/2.0,vz+h*gz2/2.0)
fy3=vy+h*gy2/2.0
gy3=accely(alpha,r,x+h*fx2/2.0,y+h*fy2/2.0,z+h*fz2/2.0, $
vx+h*gx2/2.0,vy+h*gy2/2.0,vz+h*gz2/2.0)
fz3=vz+h*gz2/2.0
gz3=accelz(alpha,r,x+h*fx2/2.0,y+h*fy2/2.0,z+h*fz2/2.0, $
vx+h*gx2/2.0,vy+h*gy2/2.0,vz+h*gz2/2.0)
fx4=vx+h*gx3
gx4=accelx(alpha,r,x+h*fx3,y+h*fy3,z+h*fz3,vx+h*gx3, $
vy+h*gy3,vz+h*gz3)
fy4=vy+h*gy3
gy4=accely(alpha,r,x+h*fx3,y+h*fy3,z+h*fz3,vx+h*gx3, $
vy+h*gy3,vz+h*gz3)
fz4=vz+h*gz3
gz4=accelz(alpha,r,x+h*fx3,y+h*fy3,z+h*fz3,vx+h*gx3, $
vy+h*gy3,vz+h*gz3)
x=x+h*(fx1+2.0*fx2+2.0*fx3+fx4)/6.0
vx=vx+h*(gx1+2.0*gx2+2.0*gx3+gx4)/6.0
y=y+h*(fy1+2.0*fy2+2.0*fy3+fy4)/6.0
vy=vy+h*(gy1+2.0*gy2+2.0*gy3+gy4)/6.0
z=z+h*(fz1+2.0*fz2+2.0*fz3+fz4)/6.0
vz=vz+h*(gz1+2.0*gz2+2.0*gz3+gz4)/6.0
r=sqrt(x*x+y*y+z*z)
x1=[x,x]
y1=[y,y]
z1=[z,z]
plots,x1,y1,z1,psym=3,color=2,/t3d
endwhile
end
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
