drawing a shaded sphere

[boccio(at)swarthmore.edu (John Boccio)]

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*rearth], $
        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