;------------------------------------------------------------ ;EXAMPLE of the driver procedure for rectifying Cassini images ;------------------------------------------------------------ program='N1458972629_1' ps=0.1 ;------------------------------------------------------------ ;using calibrated images with known scale imadir='~/CISSCAL/images/' imadir='' print,'read the image file ?' print,'y-yes cr - skip' vast='' read,vast if(vast eq 'y') then begin restore,imadir+'N1458972629_1.IMG.cal.sav' ima=rotate(writebuff,4) endif ;NA image usewa=0 nx=n_elements(ima(*,0)) ny=n_elements(ima(0,*)) ;FROM RINGS NODE: ;2004 86 5 37 48180864. 65.803 -16.319 -25.123 85.484 155.945 ;2004 86 5 38 48180542. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 39 48180219. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 40 48179897. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 41 48179575. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 42 48179253. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 43 48178930. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 44 48178608. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 45 48178286. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 46 48177963. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 47 48177641. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 48 48177319. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 49 48176997. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 50 48176674. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 51 48176352. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 52 48176030. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 53 48175708. 65.804 -16.319 -25.123 85.484 155.945 ;2004 86 5 54 48175385. 65.804 -16.319 -25.123 85.484 155.945 ;N1458972629_1.IMG ;2004-086T05:47:08.803 <--- IMAGE_TIME ;initial guesses for the rectifying process ;pivot point: use the header coordinates of the image center = planet center=1 xs0=0. ys0=0. ximage0=509.10299-1. yimage0=524.97284-1. xc_cas=0. ; assume the camera was directed to planet yc_cas=0. ; r_cas=48183.398 ;image header r_cas=48177.641 ;rings node b_cas=-16.319 ; rot_cas=0. ;BUT THIS LOOKS MUCH BETTER?? ;ARE THE PIXEL SCALES WRONG IN carrini_proje2005_f.pro ? r_cas=48300. ;making a fit with the Encke-division: ; using a +/- 5000 km wide region around the center of the Encke division ; resampling the image into 50 radial and 180 azimuthal pixels encke=5 nfit=[50,180] ;exclude from the fit: planet= [150,210], shadow [240,250] exclude=[150,210,240,270] ;---------------------------------------------------------- ;THESE INITIAL VALUES WERE USED IN THE ITERATION ;this was run, giving the following improved values ;(total squares deviation 105 --> 2.82) ; 0.11049545 508.52515 523.54588 48411.307 ; -16.310996 0.0000000 0.0000000 2.8295676 ;also, the exclude-setting were from an old file ;better one exclude=[0,5,150,210] rot_cas=0.11049545 ximage0=508.52515 yimage0=523.54588 r_cas=48411.307 b_cas=-16.310996 ;also, use abetter sampled image strip nfit=[50,180] encke=2 ;initial deviation 2.58746 was not improved ;------------------------------------------------------------------- ;BEFORE FIT: plot the Cassini image + Encke division ;------------------------------------------------------------------- print,'show initial guess ?' print,'y-yes ret-skip q-quit' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then begin loadct,3 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,rot=rot_cas,$ image=ima,pimage=[ximage0,yimage0,xc_cas,yc_cas,1.],$ grid=133.59+[-1,1]*.05 ots,'r_cas='+string(r_cas) endif ;---------------------------------------------------------------------- ;THIS MAKES THE ITERATION - AND TAKES A L O N G T I M E ;---------------------------------------------------------------------- print,'make the iteration ?' print,'y-yes ret-skip q-quit' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then begin ;what variables to vary in the fit, and on what order usefit=0 usefit=[4,5,6,0,1] p=[r_cas,b_cas,xc_cas,yc_cas,rot_cas,ximage0,yimage0] scale=1. guess=0. loadct,3 cassini_fit2005_f,ima,p,pout,scale=scale,usewa=usewa,$ exclude=exclude,usefit=usefit,guess=guess,center=center,nfit=nfit endif ;---------------------------------------------------------------------- print,'display image + grid ' print,'y-yes ret-skip q-quit' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then begin print,'image size=',nx,ny print,'give scale' read,scale nwin,xs=nx/scale,ys=ny/scale tvscl,congrid(ima,nx/scale,ny/scale) cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ inter=[80,140,14,360],usewa=usewa plots,(ximage0+xs-xs0)/scale,(yimage0+ys-ys0)/scale,$ psym=3,/dev,col=255 for i=-17,15 do begin cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,$ rot=rot_cas,$ inter=[75.,145.,700.,1.,10.*i,10*i],usewa=usewa plots,(ximage0+xs-xs0)/scale,(yimage0+ys-ys0)/scale,$ psym=3,/dev,col=255 endfor ots,program print,'dump to file?' vast='' read,vast if(vast eq 'y') then psdump,program+'_image_grid.ps' endif ;------------------------------------------------ ;slide-image view + some distances marked ? ;------------------------------------------------ print,'make the sliding image?' print,'y-yes cr-skip' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then begin print,'***making the sliding image' nwin,xs=nx,ys=ny,/pix tvscl,ima cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ inter=[70,140,7,3600],usewa=usewa plots,ximage0+(xs-xs0),yimage0+(ys-ys0),$ psym=3,/dev,col=255 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ inter=[70.,140.,7000.,1.,90.,90.],usewa=usewa plots,ximage0+(xs-xs0),yimage0+(ys-ys0),$ psym=3,/dev,col=255 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ inter=[70.,140.,700.,1.,80.,80.],usewa=usewa plots,ximage0+(xs-xs0),yimage0+(ys-ys0),$ psym=3,/dev,col=255 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ inter=[70.,140.,700.,1.,70.,70.],usewa=usewa plots,ximage0+(xs-xs0),yimage0+(ys-ys0),$ psym=3,/dev,col=255 ;camera pointing cassini_proje2005_f,xc_cas,yc_cas,xc_cas*0.,xc_cas,yc_cas,b_cas,r_cas,$ xs,ys,rot=rot_cas,usewa=usewa plots,ximage0+(xs-xs0),yimage0+(ys-ys0),$ psym=6,/dev,col=255 ima2=tvrd(0,0) slide_image,ima2,title=program+strjoin(string([p,xs0,ys0])),$ xvis=600,yvis=1.*ny/nx*600. endif ;------------------------------------------------------- ; using cassini_strip2005 for making a rectified image ;------------------------------------------------------- print,'make rectified image: rfima?' print,'y-yes ret-skip, q-end' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then begin print,'***calling cassini_strip' cassini_strip2005,ima,p,xs0,ys0,rr=[80,140],$ imaout=rfima,xima=rima,yima=fima,$ nrfii=[1201,720],/iss oplot,136.51*[1,1],[-100,100]*20,lines=2 oplot,133.43*[1,1],[-100,100]*20,lines=2 oplot,133.75*[1,1],[-100,100]*20,lines=2 oplot,84.74*[1,1],[-100,100]*20,lines=2 oplot,84.94*[1,1],[-100,100]*20,lines=2 oplot,86.37*[1,1],[-100,100]*20,lines=1 oplot,86.59*[1,1],[-100,100]*20,lines=1 save,file=program+'_rfima.save',rfima,rima,fima endif skip: ;------------------------------------------------------- ; show the rectified image + mark 4 radial regions ;------------------------------------------------------- psdirect,program+'_recti',/color,ps restore,program+'_rfima.save' nwin,xs=1000,ys=600 rfima256=rfima/max(rfima)*256 loadct,3 contour,smooth(rfima256,3),rima,fima,xr=[90,140],$ xtitle='rad',ytitle='longitude',title=program,$ lev=40+findgen(25)*4,c_col=50+findgen(25)*6,/fil contour,smooth(rfima256,5),rima,fima,xr=[90,140],/over,/noerase,$ lev=[50,55,60],c_col=255-[250,240,230] mask=rfima256 mask(800:*,*)=255 contour,smooth(mask,5),rima,fima,xr=[90,140],/over,/noerase,$ lev=[105,110,115],c_col=[250,240,230] contour,smooth(mask,5),rima,fima,xr=[90,140],/over,/noerase,$ lev=[60,65,70],c_col=255-[250,230,210] ;regions of interest: rmintab=[127.5,108.7,114.,96.5] rmaxtab=[128.5,109.7,115.,97.5] for ii=0,3 do begin oplot,rmintab(ii)*[1,1],[0,360],lines=2 oplot,rmaxtab(ii)*[1,1],[0,360],lines=2 endfor psdirect,program+'_recti',/color,ps,/stop mincurves: psdirect,program+'_mincurves',/color,ps !p.multi=[0,1,4] !p.charsize=1. nwin,xs=600,ys=800 yrange=[50,200] for ii=0,3 do begin ind=where(rima gt rmintab(ii) and rima lt rmaxtab(ii)) stab=string(rmintab(ii),'(f5.1)')+'-'+string(rmaxtab(ii),'(f5.1)') plot,fima,rfima(ind(0),*),yr=[.1,.3],ys=1,psym=3,$ title=program+stab,xtitle='longitude',ytitle='I/F' for i=0,n_elements(ind)-1 do begin oplot,fima,rfima(ind(i),*),psym=3 endfor endfor !p.multi=0 label_data,0.1,0.9,'SIN-CURVE MIN AT 110',col=160,size=1 psdirect,program+'_mincurves',/color,ps,/stop ;goto,endi eit: ;goto,skippi ;------------------------------------------------------ minima: restore,program+'_rfima.save' rfima256=rfima/max(rfima)*256. psdirect,program+'_minima',/color,ps nwin,xs=1000,ys=600 loadct,3 contour,smooth(rfima256,3),rima,fima,xr=[90,140],$ xtitle='rad',ytitle='longitude',title=program,$ lev=50+findgen(25)*4,c_col=50+findgen(25)*6,/fil n800=200 n1200=1149 calc=0 print,'calculate minima (y, ret-no, q-quit)' vast='' read,vast if(vast eq 'q') then goto,endi if(vast eq 'y') then calc=1 if(calc eq 1) then begin mintab2=rima*0. amp20=mintab2*0. amp10=amp20 amp30=amp20 endif if(calc eq 0) then begin restore,program+'_rfima_min.save' endif fiirange=where(fima gt 55 and fima lt 140) for i=n800,n_elements(rima)-1 do begin apu=rfima(i,fiirange) ;if(rima(i) gt 122 and rima(i) le 137.) then begin if(rima(i) gt -122 and rima(i) le 100+137.) then begin if(calc eq 1) then begin mc_fit_minimum,fima(fiirange),reform(apu),$ lmin=lmin,amp_20=amp_20,amp_30=amp_30,amp_10=amp_10 mintab2(i)=lmin amp10(i)=amp_10-1. amp20(i)=amp_20-1. amp30(i)=amp_30-1. endif if(amp20(i) gt 0) then begin plots,rima(i),mintab2(i),psym=1,col=220 endif endif endfor skippi: if(calc eq 1) then begin save,file=program+'_rfima_min.save',$ rfima,rima,fima,amp10,amp20,amp30,mintab2 endif ind=where(RIMA LT 117 OR (rima gt 122 and rima lt 137)) oplot,rima(ind),mintab2(ind),psym=1,col=20 oplot,rima,rima*0+90,lines=2 oplot,rima,rima*0+100,lines=2 psdirect,program+'_minima',/color,ps,/stop psdirect,program+'_min_vs_tau',/color,ps restore,'/home/heikki/MCCODE/RESULTS/HST/tau_pps_smoothed.sav',/v rtab=r_vals/1000. tautab=tau_pps_smoothed y2=spl_init(rtab,tautab) tauima=spl_interp(rtab,tautab,y2,rima) nwin !p.multi=[0,2,1] tek_color plot,tauima,mintab2,/nod,xtitle='TAU',ytitle='minimum longitude',title=program,xr=[0,2.5],xs=1,yr=[50,150],ys=1 ind=where(rima gt 93 and rima lt 98) oplot,tauima(ind),mintab2(ind),psym=6,sym=.5,col=2 ind=where(rima gt 98 and rima lt 105) oplot,tauima(ind),mintab2(ind),psym=6,sym=.5,col=3 ind=where(rima gt 105 and rima lt 115) oplot,tauima(ind),mintab2(ind),psym=4,sym=.5,col=4 ind=where(rima gt 123 and rima lt 131) oplot,tauima(ind),mintab2(ind),psym=1,col=5 label_data,0.01,0.9,['93-98','98-105','105-115','123-131'],col=[2,3,4,5],psym=[6,6,4,1],syms=[.5,.5,.5,1] n90=n_elements(rfima(0,*))/4 iftab2=rfima(*,n90) plot,tauima,iftab2,/nod,xtitle='PPS TAU',ytitle='ADU',title=program,xr=[0,2.5],xs=1,yr=[50,150]*0,ys=1 ind=where(rima gt 93 and rima lt 98) oplot,tauima(ind),iftab2(ind),psym=6,sym=.5,col=2 ind=where(rima gt 98 and rima lt 105) oplot,tauima(ind),iftab2(ind),psym=6,sym=.5,col=3 ind=where(rima gt 105 and rima lt 115) oplot,tauima(ind),iftab2(ind),psym=4,sym=.5,col=4 ind=where(rima gt 123 and rima lt 131) oplot,tauima(ind),iftab2(ind),psym=1,col=5 !p.multi=0 psdirect,program+'_min_vs_tau',/color,ps,/stop goto,endi ind=where((rima gt 122 and rima lt 137) AND AMP20 GT 0) oplot,rima(ind),amp20(ind)*200,psym=1,col=220 ;,NSUM=5 oplot,rima(ind),amp30(ind)*200,psym=3,col=220 ;,NSUM=5 oplot,rima(ind),amp10(ind)*200,psym=3,col=220 ;,NSUM=5 ind=where(RIMA LT 117 AND AMP20 GT 0 AND MINTAB2 LE 100) oplot,rima(ind),amp20(ind)*200,psym=1,col=220 ;,NSUM=5 oplot,[120,140],[40,40],col=255 oplot,[1,1]*122,[0,40-1],col=255 oplot,[0,-0.2]+122,[20,20],col=255 xyouts,120.5,10,ori=90,'A!b10!n A!b20!n A!b30!n',col=255 xyouts,121.2,19,'0.1',col=255 endi: end