pro cassini_fit2005_f,image,pin,pout,keeler=keeler,cring=cring,encke0=encke0,$ usewa=usewa,iplot0=iplot0,r00=r00,scale0=scale0,nfit=nfit,$ usefit0=usefit0,exclude0=exclude0,guess=guess,$ center0=center0,$ rfit0=rfit0, rmin_fit0=rmin_fit0, rmax_fit0=rmax_fit0,$ fiimin_fit0=fiimin_fit0, fiimax_fit0=fiimax_fit0 common cassini_model_com1,ima,$ rmin_fit,rmax_fit,fiimin_fit,fiimax_fit,nr_fit,nfii_fit,rfit common cassini_model_com2,iplot common cassini_model_com3,r0,fii0,xs0,ys0 common cassini_model_com4,ite_counter common cassini_model_com5,exclude, pfix, usefit common cassini_model_com6,center common cassini_model_com7,usewa0 ;related to the pivot point, see below center=0 if(keyword_set(center0)) then center=1 ;exclude somelongitudes from the fit? exclude=0 if(keyword_set(exclude0)) then exclude=exclude0 ;----------------------------------------------------------------- ;initial paramaters for the iteration ;----------------------------------------------------------------- pfix=pin r_cas=pin(0) b_cas=pin(1) xc_cas=pin(2) yc_cas=pin(3) rot_cas=pin(4) ximage0=pin(5) yimage0=pin(6) ;specify the paramaters to be fitted = indexes in usefit ;others are kept at their initial values usefit=[0,1,2,3,4,5,6] if(keyword_set(usefit0)) then usefit=usefit0 p=dblarr(n_elements(usefit)) p_ide='' for i=0,n_elements(usefit)-1 do begin if(usefit(i) eq 0) then p(i)=r_cas if(usefit(i) eq 1) then p(i)=b_cas if(usefit(i) eq 2) then p(i)=xc_cas if(usefit(i) eq 3) then p(i)=yc_cas if(usefit(i) eq 4) then p(i)=rot_cas if(usefit(i) eq 5) then p(i)=ximage0 if(usefit(i) eq 6) then p(i)=yimage0 if(usefit(i) eq 0) then p_ide=p_ide+' r_cas' if(usefit(i) eq 1) then p_ide=p_ide+' b_cas' if(usefit(i) eq 2) then p_ide=p_ide+' xc_cas' if(usefit(i) eq 3) then p_ide=p_ide+' yc_cas' if(usefit(i) eq 4) then p_ide=p_ide+' rot_cas' if(usefit(i) eq 5) then p_ide=p_ide+' ximage0' if(usefit(i) eq 6) then p_ide=p_ide+' yimage0' endfor p_ide=p_ide+' xs0 ys0 DEV' ;counter for calls to cassini_model ite_counter=0 ;default pixels size = Cassini NA pixel usewa0=0 if(keyword_set(usewa)) then usewa0=1 ;cassini image ima=image nx=n_elements(ima(*,0)) ny=n_elements(ima(0,*)) nxc0=0.5*double(nx-1.d0) nyc0=0.5*double(ny-1.d0) ;scale reduction factor for plotting ima scale=5 if(keyword_set(scale0)) then scale=scale0 ;----------------------------------------------------------------- ;defaults for the accuracy of the resampled image region fitted ;----------------------------------------------------------------- nr_fit=50 nfii_fit=180 if(keyword_set(nfit)) then begin nr_fit=nfit(0) nfii_fit=nfit(1) endif ;-------------------------------------------------------------------- ;Define: ;Encke division center: rfit ;region around Encke division: rmin rmax ;longitude range of fit: fiimin,fiimax ;-------------------------------------------------------------------- encke=1.d0 if(keyword_set(encke0)) then encke=encke0 rfit=133.59d0 rmin_fit=133.d0-encke rmax_fit=134.d0+encke fiimin_fit=0.d0 fiimax_fit=360.d0 ;add Keelere gap to the fit? if(keyword_set(keeler)) then begin rfit =[rfit, 136.509d0] rmin_fit =[rmin_fit, 136.200d0] rmax_fit =[rmax_fit, 136.650d0] fiimin_fit =[fiimin_fit, 70.d0 ] fiimax_fit =[fiimax_fit, 110.d0 ] endif ;take from c-ring: max -neg? if(keyword_set(cring)) then begin rfit =[rfit, -86.48d0] rmin_fit =[rmin_fit, 85.50d0] rmax_fit =[rmax_fit, 87.50d0] fiimin_fit =[fiimin_fit, fiimin_fit] fiimax_fit =[fiimax_fit, fiimax_fit] endif if(keyword_set(rfit0)) then rfit=rfit0 if(keyword_set(rmin_fit0)) then rmin_fit=rmin_fit0 if(keyword_set(rmax_fit0)) then rmax_fit=rmax_fit0 if(keyword_set(fiimin_fit0)) then fiimin_fit=fiimin_fit0 if(keyword_set(fiimax_fit0)) then fiimax_fit=fiimax_fit0 ;-------------------------------------------------------------- ;Pivot point = viewing point xs0,ys0 corresponding to ; image pixel ximage0, yimage0 ;/center -> use planet center as a pivotpoint xs0=ys0=0 ;no /center -> use the ansa point at some distance r0 ; default= ansa point of Encke division inner edge r0=133.425 if(keyword_set(r00)) then r0=r00 if(center eq 0) then begin ;-------------------------------------------------------------------- ;find the apparant ansa point at a given distance r0 ;-------------------------------------------------------------------- nfii=360 & fiimin=30. & fiimax=150. fiit=(dindgen(nfii*10l+1l)/(nfii*10l)*(fiimax-fiimin)+fiimin)*!dpi/180.d0 x0=[0,r0*cos(fiit)] y0=[0,r0*sin(fiit)] z0=x0*0.d0 cassini_proje2005_f,x0,y0,z0,xc_cas,yc_cas,b_cas,r_cas,xs00,ys00,$ rot=rot_cas,usewa=usewa ;find the apparant ansa - maximum distance from the planet xs_plan=xs00(0) ys_plan=ys00(0) xs00=xs00(1:*) ys00=ys00(1:*) rtemp=sqrt((xs00-xs_plan)^2+(ys00-ys_plan)^2) ind=where(rtemp eq max(rtemp)) ind=ind(0) xs0=xs00(ind) ys0=ys00(ind) print,'------------------------------' print,'XS0,YS0',xs0,ys0 print,'------------------------------' endif if(center eq 1) then begin ;-------------------------------------------------------------------- ;center image at the planet (instead of the pivot point) ;-------------------------------------------------------------------- xs0=0.d0 ys0=0.d0 endif ;-------------------------------------------------------------------- ; MISCELLENOUS CHECKS ;-------------------------------------------------------------------- ;-------------------------------------------------------- ;plot the ring grid corresponding to initial paramaters cassini_proje2005_f,x0,y0,z0,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ grid=[70,137.6,32,360],mag=5,/plot,usewa=usewa oplot,[-1,-1,1,1,-1]*512.,[-1,1,1,-1,-1]*512 vast='' print,'hit cr to continue...' read,vast wdelete ;-------------------------------------------------------- ;plot the image + initial grid nwin,xs=nx/scale,ys=ny/scale tvscl,congrid(ima,nx/scale,ny/scale) nfii=65 fiimax=fiimax_fit(0) fiimin=fiimin_fit(0) fii=(dindgen(nfii+1)/nfii*(fiimax-fiimin)+fiimin)*!dpi/180.d0 sinf=sin(fii) cosf=cos(fii) ;Encke gap - center x=133.59*cosf y=133.59*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=6,/dev,col=120,sym=1 ;Encke gap - inner x=133.43*cosf y=133.43*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=255 ;Encke gap - outer x=133.75*cosf y=133.75*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=255 ;Keeler x=136.51*cosf y=136.51*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=255 ;C ring ringlet - inner x=86.37*cosf y=86.37*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=255 ;C ring ringlet - outer x=86.59*cosf y=86.59*sinf z=x*0.d0 cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,rot=rot_cas,$ usewa=usewa plots,(ximage0+(xs-xs0))/scale,(yimage0+(ys-ys0))/scale,$ psym=3,/dev,col=255 ;Pivot point 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))/scale,(yimage0+(ys-ys0))/scale,$ psym=6,/dev,col=255 ots,'Initial guess for the grid' vast='' print,'hit cr to continue...' read,vast ;---------------------------------------- ; rectified image cassini_strip2005,ima,[r_cas,b_cas,xc_cas,yc_cas,rot_cas,ximage0,yimage0],$ xs0,ys0,rr=[80.,140.],nrfii=[600,100],$ ima=rfima,xima=rima,yima=fima,/iss,usewa=usewa 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 if(keyword_set(guess)) then begin new_guess: print,'current values:' print,pfix print,'ide=r,b,xc,yc,rot,x0,y0 - > modify this' print,'ide='' -> proceed to fit' ide='' read,ide,val if(ide eq 'r') then r_cas=val if(ide eq 'b') then b_cas=val if(ide eq 'xc') then xc_cas=val if(ide eq 'yc') then yc_cas=val if(ide eq 'rot') then rot_cas=val if(ide eq 'x0') then ximage0=val if(ide eq 'y0') then yimage0=val if(ide eq '') then goto,end_guess pfix=[r_cas,b_cas,xc_cas,yc_cas,rot_cas,ximage0,yimage0] for i=0,n_elements(usefit)-1 do begin if(usefit(i) eq 0) then p(i)=r_cas if(usefit(i) eq 1) then p(i)=b_cas if(usefit(i) eq 2) then p(i)=xc_cas if(usefit(i) eq 3) then p(i)=yc_cas if(usefit(i) eq 4) then p(i)=rot_cas if(usefit(i) eq 5) then p(i)=ximage0 if(usefit(i) eq 6) then p(i)=yimage0 endfor cassini_strip2005,ima,[r_cas,b_cas,xc_cas,yc_cas,rot_cas,ximage0,yimage0],$ xs0,ys0,rr=[80.,140.],nrfii=[600,100],$ ima=rfima,xima=rima,yima=fima,/iss,usewa=usewa 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 goto, new_guess endif end_guess: ;------------------------------------------------------------ ;start the actual fit with Powell method ;------------------------------------------------------------ np=n_elements(usefit) xi=dblarr(np,np) for i=0,np-1 do begin for j=0,np-1 do begin if(i eq j) then xi(i,j)=1.0d0 endfor endfor itmax=20 ftol=1d-5 nwin !p.position=[0.1,0.1,0.7,0.95] iplot=1 junk=cassini_model2005(p) vast='' print,'hit cr to continue...' read,vast itmax=2000 iplot=10 print,p_ide powell,p,xi,ftol,fmin,'cassini_model2005',itmax=itmax,iter=iter,/double print,p !p.position=0 pout=p endi: end