;---------------------------------------------------------------- function cassini_model2005,p,plot=plot ;---------------------------------------------------------------- If(n_params() le 0) then begin print,'-------------------------------------------------------------------' print,'function cassini_model2005,p,plot=plot' print,'-------------------------------------------------------------------' print,'uses cassini_proje2005_f to produce a grid for a ring region' print,' [RMIN, RMAX, FIIMIN, FIIMAX] (in xxx_fit common variables)' print,'returns the deviation from the center of the Encke gap' print,' or some other circular feature defined via variable RFIT (common)' print,' --> used as a function to be minimized in cassini2005_fit' print,'-------------------------------------------------------------------' print,' Model parameters:' print,' p=[R_CAS, B_CAS , XC_CAS, YC_CAS, ROT_CAS, XIMAGE0, YIMAGE0]' print,'' print,' R_CAS = distance of camera from Saturn [1000 km]' print,' B_CAS = elevation angle (degrees)' print,' XC_CAS YC_CAS = ring point to which the camera is directed' print,' x-axis toward sub-observer point [1000 km]' print,' ROT_CAS = CCW rotation of camera field [degrees]' print,' wrt to ori with local North pointing up' print,' XIMAGE0,YIMAGEI0 = pixel location in image IMA (common)' print,' corresponding to a ring view location XS0,YS0' print,' XS0, YS0 can be defined directly or via a ring location' print,' R0,FII0 (i.e. inner edge of Encke at ansa) (common)' print,'-------------------------------------------------------------------' print,'' print,' procedure cassini2005_proje_f calculates:' print,' x,y,z -> xs, ys arbitrary x,y,z with given R,B,XC,YC,ROT' print,' x0,y0,z0 -> xs0,ys0 for r0,fii0 point -"-' print,' print,' define the XIMAGE,YIMAGE grid on top of image IMA:' print,' XIMAGE - XIMAGE0 = (XS-YS0)*scale' print,' YIMAGE - YIMAGE0 = (YS-YS0)*scale' print,' default scale is that of Cassini NA camera' print,'-------------------------------------------------------------------' return,'' endif 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 ite_counter=ite_counter+1 usewa=usewa0 ;----------------------------------------------------------------- ;usefit contains indices of variables included to the fit ;others are kept fixed to their initial values r_cas=pfix(0) b_cas=pfix(1) xc_cas=pfix(2) yc_cas=pfix(3) rot_cas=pfix(4) ximage0=pfix(5) yimage0=pfix(6) for i=0,n_elements(usefit)-1 do begin if(usefit(i) eq 0) then r_cas=p(i) if(usefit(i) eq 1) then b_cas=p(i) if(usefit(i) eq 2) then xc_cas=p(i) if(usefit(i) eq 3) then yc_cas=p(i) if(usefit(i) eq 4) then rot_cas=p(i) if(usefit(i) eq 5) then ximage0=p(i) if(usefit(i) eq 6) then yimage0=p(i) endfor nx=n_elements(ima(*,0)) ny=n_elements(ima(0,*)) nxc=0.5d0*(nx-1.d0) nyc=0.5d0*(ny-1.d0) nfii=nfii_fit fiimin=fiimin_fit(0) fiimax=fiimax_fit(0) ;---------------------------------------------------------------- ;What is the pivot point xs0,ys0 corresponding to XIMAGE0,YIMAGE0 ;---------------------------------------------------------------- ;determined by center=1 or 0 ;--------------------------------------------- ;center image at the planet ? if(center eq 1) then begin xs0=0.d0 ys0=0.d0 endif ;--------------------------------------------- ;the apparant ansa point at a given distance r0 ;(so fii0 is ignored) if(center ne 1) then begin fiit=(dindgen(nfii*10l+1l)/(nfii*10l)*(fiimax-fiimin)+fiimin)*!dpi/180.d0 x0=[0.d0,r0*cos(fiit)] y0=[0.d0,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 ;find ansa? 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) endif ;-------------------------------------------------------------- ;make an image of the fitting region ;from which the deviattion is calculated ;-------------------------------------------------------------- ;either 1 or 2 fitted regions -> plot both if needed nfit=n_elements(rfit) if(nfit eq 1) then begin pos=fltarr(1,4) pos(0,*)=[0.1,0.1,0.95,0.9] endif if(nfit eq 2) then begin pos=fltarr(2,4) pos(0,*)=[0.05,0.1,0.45,0.9] pos(1,*)=[0.50,0.1,0.95,0.9] endif ;---------------------------------------------- ;loop over fitting regions devi=0. for i=0,nfit-1 do begin ;set inter to the paramaters of the fitting region r1=rmin_fit(i) r2=rmax_fit(i) fiimin=fiimin_fit(i) fiimax=fiimax_fit(i) nr=nr_fit inter=[r1,r2,nr,nfii,fiimin-10.d0,fiimax+10.d0] fiid=(dindgen(nfii+1)/nfii*(fiimax-fiimin+20.d0)+fiimin-10.d0) fii=fiid*!dpi/180.d0 rad=(dindgen(nr+1)/nr*(r2-r1)+r1) cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,$ inter=inter,rot=rot_cas,usewa=usewa xint=ximage0+(xs-xs0) yint=yimage0+(ys-ys0) ;interpolate image values corresponding to the fitting region mask=fltarr(nfii+1,nr+1) image_bilin,ima,xint,yint,adut mask(0:nfii,0:nr)=adut ;at each azimuth, find the minimum brightness ;within the azimuthal limits of the fitting region ind_fii_use=where(fiid ge fiimin and fiid le fiimax,count) ;exclude the planet shining through ring, shadow? ;exclude-keyword, exclude maximum of 3 longitude intervals: ;exclude=[fiimin1,fiimax1, fiimin2,fiimax2, fiimin3,fiimax3] if(n_elements(exclude) ge 2) then begin ind_not_exclude= where(fiid le exclude(0) or fiid ge exclude(1)) ind_fii_use=common_index(ind_fii_use,ind_not_exclude) if(n_elements(exclude) ge 4) then begin ind_not_exclude= where(fiid le exclude(2) or fiid ge exclude(3)) ind_fii_use=common_index(ind_fii_use,ind_not_exclude) endif if(n_elements(exclude) ge 6) then begin ind_not_exclude= where(fiid le exclude(4) or fiid ge exclude(5)) ind_fii_use=common_index(ind_fii_use,ind_not_exclude) endif endif count=n_elements(ind_fii_use) rminloc=fltarr(count) ;bad = number of azimuths where brightnes minimum falls on ;the radial boundaries bad=0. ;rfit = positive indicates finding a minimum (for a division) IF(RFIT(I) GT 0) THEN BEGIN for ix=0,count-1 do begin iix=ind_fii_use(ix) indmin=where(mask(iix,*) eq min(mask(iix,*))) rminloc(ix)=rad(indmin(0)) if(indmin(0) eq 0 or indmin(0) ge nr) then bad=bad+1. endfor ENDIF ;rfit = negative indicates finding a maximum (for a ringlet) IF(RFIT(I) LT 0) THEN BEGIN for ix=0,count-1 do begin iix=ind_fii_use(ix) indmin=where(mask(iix,*) eq max(mask(iix,*))) rminloc(ix)=rad(indmin(0)) if(indmin(0) eq 0 or indmin(0) ge nr) then bad=bad+1. endfor ENDIF ;deviation from the center location: ;total of squares, taking into account bad's devi=devi+total((rminloc-abs(rfit(i)))^2)*(1.d0+bad/float(i+1)) ;-------------------------------------------------------- ;plot? if(iplot ne 0) then begin !p.position=pos(i,*) tvplot2,mask,xx=fii*!radeg,yy=rad,/nonew,noeras=i oplot,fii*!radeg,fii*0+abs(rfit(i)),lines=2 oplot,fii(ind_fii_use)*!radeg,rminloc,psym=6,thick=3 oplot,fiimin*[1,1],[0,1000],lines=1 oplot,fiimax*[1,1],[0,1000],lines=1 if(n_elements(exclude) ge 2) then begin oplot,exclude(0)*[1,1],[0,1000],lines=2,thick=2 oplot,exclude(1)*[1,1],[0,1000],lines=2,thick=2 endif if(n_elements(exclude) ge 4) then begin oplot,exclude(2)*[1,1],[0,1000],lines=3 oplot,exclude(3)*[1,1],[0,1000],lines=3 endif if(n_elements(exclude) ge 6) then begin oplot,exclude(4)*[1,1],[0,1000],lines=3 oplot,exclude(5)*[1,1],[0,1000],lines=3 endif !p.position=0 endif endfor if(iplot ge 1) then begin ff1='(f8.3)' ff0='(f8.2)' ff2='(f6.2)' ff5='(f9.6)' title='DIST='+string(r_cas,ff1)+' B='+string(b_cas,ff2)+$ ' XC,YC='+string(xc_cas,ff0)+string(yc_cas,ff0)+$ ' rot='+string(rot_cas,ff0)+'!c!c'+$ ' XI0,YI0='+string(ximage0,ff0)+string(yimage0,ff0)+$ ' dev='+string(devi,ff5)+' bad='+string(bad,ff2) ots,title,size=1,/arp wait,.3 endif if(iplot gt 1) then begin if(ite_counter ge iplot) then begin erase 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,/nonew,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 vast='' ; read,vast wait,2 ite_counter=0 endif endif title=strjoin([strtrim(p),string(xs0),string(ys0),string(devi)]) print,title return,devi end