pro cassini_proje2005_f,x0,y0,z0,xc0,yc0,b_cas0,r_cas0,rot0=rot0,xs,ys,$ plot=plot,test=test,grid=grid,autorot0=autorot0,$ pitch=pitch,usewa=usewa,$ nonew=nonew,zc0=zc0,xorb=xorb,yorb=yorb,zorb=zorb,$ add=add,inter=inter,$ magfac0=magfac0,image=image,pimage=pimage if(n_params() le 0) then begin print,'-----------------------------------------------------------' print,'cassini_proje2005_f,x,y,z,xc,yc,b_cas,r_cas,xs,ys,rot0=rot0' print,'-----------------------------------------------------------' print,'transforms 3D points (x,y,z) to Cassini image plane (xs,ys)' print,'checked 27.10.2005 HS' PRINT,' ' print,'INPUT:' print,' x,y,z table of 3D points viewed from Cassini' print,' xc,yc = intersection of viewing direction in the ring plane' print,' in complete images this is usually =0,0?' print,' [coordinate system: x=sub-observer, z = North] print,' b_cas = elevation (in degrees) of Cassini wrt Saturn center' print,' r_cas = distance of Cassini from Saturn center [unit= 1000 km]' print,'OUTPUT:' print,' projected points xs,ys in NA camera pixels' print,'KEYWORDS:' print,' useWA -> calculate xs,ys,assuming WA-camera' print,' rot = apply ccw rotation to the final image' print,' xorb,yorb,zorb -> use these instead r_cas,b_cas' print,' zs0 -> use this for zs (def=0.)' print,'' print,' /plot' print,' magfac=value -> expand the plot region (default= Cassini view)' print,' /test -> add local ring plane normal at view center' print,' /autorot -> rotate so that local normal at 0,0 is perpendicular' print,' /grid -> create standard grid of values x,y,z' print,' or grid=[rmin,rmax,nr,nfii] print,' pitch=value -> add wakes to the grid' print,' /add -> add 0,0,0 to grid print,' /inter-> create standard grid used in interpolation' print,' or inter=[rmin,rmax,nr,nfii,fiimin,fiimax]' print,'' print,' image=image and pimage=[ximage0,yimage0,xs0,ys0,scale] given' print,' --> display the points [xs,ys] on top of the Cassini image' print,' ximage0,yimage0 is the pixel location in the image' print,' corresponding to a specific xs0,ys0 (usually planet center)' print,' scale = factor to scale the plot' print,'-----------------------------------------------------------' print,'examples' print,';NA view from 48 million kilometers' print,'cassini_proje2005_f,x,y,z,0,0,20.,48000.,/grid,/plot' print,';100*NA view from 0.48 million kilometers (illustrate perspective)' print,'cassini_proje2005_f,x,y,z,0,0,20.,480.,/grid,/plot,mag=100' print,';same toward ansa print,'cassini_proje2005_f,x,y,z,0,100,20.,480.,/grid,/plot,mag=100,/test' print,'cassini_proje2005_f,x,y,z,0,100,20.,480.,/gr,/pl,ma=100,/te,/auto' print,';movie:' print,'xorb=500.-findgen(200)*5 & yorb=xorb*0.+100. & zorb=xorb*0.+100' print,'for i=0,199 do begin & cassini_proje2005_f,x,y,z,0,0,1.,1.,xo=xorb(i),yo=yorb(i),zo=zorb(i),nonew=i,/grid,/plot,mag=1000.,/test & endfor' return endif rad2deg=180.d0/!dpi deg2rad=!dpi/180.d0 ;-------------------------------------------------------------------------- ;use the x,y,z points given via call? ;-------------------------------------------------------------------------- ;convert to double if(not keyword_set(grid) and not keyword_set(inter)) then begin x=1.d0*x0 y=1.d0*y0 z=1.d0*z0 endif ;-------------------------------------------------------------------------- ;x,y,z defined via inter? if(keyword_set(inter)) then begin ;-------------------------------------------------------------------------- r1=60.d0 r2=200.d0 nr=121 nfii=3600 fiimin=0.d0 fiimax=360.d0 if(n_elements(inter) ge 2) then begin r1=inter(0)*1.d0 r2=inter(1)*1.d0 endif if(n_elements(inter) ge 3) then nr=inter(2) if(n_elements(inter) ge 4) then nfii=inter(3) if(n_elements(inter) ge 5) then fiimin=inter(4)*1.d0 if(n_elements(inter) ge 6) then fiimax=inter(5)*1.d0 fii=((dindgen(nfii+1))/nfii*(fiimax-fiimin)+fiimin)*deg2rad sinff=sin(fii) cosff=cos(fii) x=r1*cosff y=r1*sinff z=x*0.d0 if(nr gt 0) then begin for ir=1,nr do begin r=r1+ir*1.d0/nr*(r2-r1) x=[x,r*cosff] y=[y,r*sinff] z=[z,sinff*0.d0] endfor endif x0=x y0=y z0=z endif ;-------------------------------------------------------------------------- ;x,y,z defined via grid? if(keyword_set(grid) and not keyword_set(inter)) then begin ;-------------------------------------------------------------------------- r1=70. r2=145. nr=15 ;5000 km intervals nfii=36*5 ;circles drawn with 2 degree resolution nspoke=36 ;spoke every 10 degrees if(n_elements(grid) ge 2) then begin r1=grid(0) r2=grid(1) endif if(n_elements(grid) ge 3) then nr=grid(2) if(n_elements(grid) ge 4) then nfii=grid(3) nline=nr*5 fii=dindgen(nfii+1)/nfii*!dpi*2.d0 sinff=sin(fii) cosff=cos(fii) line=r1+dindgen(nline+1)/nline*(r2-r1) line_ansa=r1+dindgen(nline*1.1+1)/nline*(r2-r1) ;add center ? x=0.d0 y=0.d0 z=0.d0 if(keyword_set(add)) then begin x=[x,1.d0*x0] y=[y,1.d0*y0] z=[z,1.d0*z0] endif ;create nr rings for ir=0,nr-1 do begin r=r1+ir*1.d0/nr*(r2-r1) x=[x,r*cosff] y=[y,r*sinff] z=[z,sinff*0.d0] endfor ;create nspoke spokes for is=0,nspoke-1 do begin angle=is*1.d0/nspoke*2.d*!dpi line_use=line if(abs(angle-!dpi/2.) le 0.001) then line_use=line_ansa if(abs(angle-!dpi*3./2.) le 0.001) then line_use=line_ansa x=[x,line_use*cos(angle)] y=[y,line_use*sin(angle)] z=[z,line_use*0.] endfor ;add wakes? if(keyword_set(pitch)) then begin line=dindgen(101)/100.*10. for is=0,nspoke-1 do begin angle=is*1.d0/nspoke*2.d*!dpi angle2=angle-(90.-pitch)*deg2rad xt=cos(angle)*r2+line*cos(angle2) yt=sin(angle)*r2+line*sin(angle2) x=[x,xt] y=[y,yt] z=[z,xt*0] endfor endif endif ;-------------------------------------------------------------------------- xc=1.d0*xc0 yc=1.d0*yc0 zc=0.d0 if(keyword_set(zc0)) then zc=1.*zc0 b_cas=1.d0*b_cas0 r_cas=1.d0*r_cas0 ;add local normal to x,y,z, points? if(keyword_set(test)) then begin line=dindgen(101)/100.*100.*test x=[x,xc+line*0] y=[y,yc+line*0] z=[z,zc+line] x=[x,line*0] y=[y,line*0] z=[z,line*.5] endif ;-------------------------------------------------------------------------- ;HERE IS THE ACTUAL PROJECTION ;-------------------------------------------------------------------------- ;------------------------------- ;viewing direction defined by ;------------------------------- ;xc,yc point in the ring plane which maps ;to the zenter of the image ;------------------------------- ;spacecraft location defined by ;------------------------------- ; distance r_cas ; elevation B_cas ;used coordinate system: x-axis points to the projection of spacecraft ;onto the ring plane, z-axis toward pole of the rings x_cas= r_cas*cos(B_cas*deg2rad) y_cas= 0.d0 z_cas= r_cas*sin(B_cas*deg2rad) ;or given directly by if(keyword_set(xorb)) then x_cas=xorb*1.d0 if(keyword_set(yorb)) then y_cas=yorb*1.d0 if(keyword_set(zorb)) then z_cas=zorb*1.d0 ;-------------------------------- ;move into a coordinate system centered at xc,yc,zc ;calculate spacecraft location in this system ;x_cas-xc = delta *cos(B1) *cos(theta1) ;y_cas-yc = delta *cos(B1) *sin(theta1) ;z_cas-zc = delta *sin(B1) x1=x_cas-xc y1=y_cas-yc z1=z_cas-zc delta=sqrt(x1^2+y1^2+z1^2) b1=asin(z1/delta) cosb1=cos(b1) theta1=atan(y1/delta/cosb1,x1/delta/cosb1) ;ARE THESE PIXEL SIZES REALLY CORRECT- ;I HAVE FORGOT FROM WHERE THEY WERE TAKEN! napix=512.d0/3.054326190990d-3 wapix=512.d0/3.03687290d-2 ;fovinfo.txt ;CASSINI_ISS_NAC RECTANGULAR +3.054326190990 +3.054326190990 ;CASSINI_ISS_NAC_RAD CIRCULAR +1.57079633E+03 +1.57079633E+03 ;CASSINI_ISS_WAC RECTANGULAR +3.03687290E+01 +3.03687290E+01 ;CASSINI_ISS_WAC_RAD CIRCULAR +1.57079633E+03 +1.57079633E+03 dscale=napix if(keyword_set(usewa)) then dscale=wapix sinf=sin(theta1) cosf=cos(theta1) sina=sin(!dpi/2.d0-b1) cosa=cos(!dpi/2.d0-b1) dorigo=delta ;------------------------------------------- xa=x-xc ya=y-yc za=z-zc xe=-xa*sinf+ya*cosf ye=-(xa*cosf+ya*sinf)*cosa+za*sina ze= (xa*cosf+ya*sinf)*sina+za*cosa-dorigo ind=where(ze lt 0) ; only points behind the screen xe=xe(ind) ye=ye(ind) ze=ze(ind) dist=ze zs=ze+dorigo xs=-dscale*(xe/ze) ys=-dscale*(ye/ze) ;------------------------------------------- ;rotate xs,ys coordinates by the angle rot ccw ? if(keyword_set(rot0)) then begin xsmui=xs ysmui=ys sinr=sin(rot0*deg2rad) cosr=cos(rot0*deg2rad) xs= xsmui*cosr-ysmui*sinr ys= xsmui*sinr+ysmui*cosr endif ;define local normal at 0,0,0 ;------------------------------------------------------------------ ;autorot? ;search for a rotation which makes the ring normal at 0,0, vertical if(keyword_set(autorot0)) then begin xauto=[0,0]*1.d0 yauto=[0,0]*1.d0 zauto=[0,50.]*1.d0 xa=xauto-xc ya=yauto-yc za=zauto-zc xe=-xa*sinf+ya*cosf ye=-(xa*cosf+ya*sinf)*cosa+za*sina ze= (xa*cosf+ya*sinf)*sina+za*cosa-dorigo dist=ze zs_auto=ze+dorigo xs_auto=-dscale*(xe/ze) ys_auto=-dscale*(ye/ze) autorot=90.d0-atan(ys_auto(1)-ys_auto(0),xs_auto(1)-xs_auto(0))*rad2deg ;apply this rotation xsmui=xs ysmui=ys sinr=sin(autorot*deg2rad) cosr=cos(autorot*deg2rad) xs= xsmui*cosr-ysmui*sinr ys= xsmui*sinr+ysmui*cosr endif ;---------------------------------------------------------------------- ;PLOT ;---------------------------------------------------------------------- ff1='(f8.1)' ff0='(f6.0)' ff2='(f6.2)' i3='(i3)' if(keyword_set(plot)) then begin if(not keyword_set(nonew)) then nwin title='DIST='+string(r_cas,ff1)+' *1000 km B='+string(b_cas,ff2)+$ ' XC,YC='+string(xc,ff0)+string(yc,ff0) xtitle='xs (NA pixels)' ytitle='ys (NA pixels)' if(keyword_set(usewa)) then begin xtitle='xs (WA pixels)' ytitle='ys (WA pixels)' endif magfac=1.d0 if(keyword_set(magfac0)) then begin magfac=magfac*magfac0 title=title+'MAG='+string(magfac,'(f8.3)') endif plot,xs,ys,psym=3,/iso,xs=1,ys=1,xtitle=xtitle,ytitle=ytitle,$ xr=[-512,512]*magfac,yr=[-512,512]*magfac,title=title endif ;---------------------------------------------------------------------- ;PLOT ON IMAGE ;---------------------------------------------------------------------- if(keyword_set(image)) then begin nx=n_elements(image(*,0)) ny=n_elements(image(0,*)) ximage0=pimage(0) yimage0=pimage(1) xs0=pimage(2) ys0=pimage(3) scale=pimage(4) nwin,xs=nx/scale,ys=ny/scale tvscl,congrid(image,nx/scale,ny/scale) plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=120,sym=1 endif end