pro cassini_strip2005,ima,p,xs0,ys0,rr=rr,imaout=imaout,xima=xima,yima=yima,$
nrfii=nrfii,nonew=nonew,iss=iss,usewa=usewa,ff=ff,encke=encke


;this program rectifies the image, using the projection paramaters
; given in p and xs0,ys0


  r_cas=p(0)
  b_cas=p(1)
  xc_cas=p(2)
  yc_cas=p(3)
  rot_cas=p(4)
  ximage0=p(5)
  yimage0=p(6)

;-----------------------------------------------------------
;sample into nr*nfii image
;from r1 to r2   anf fiimin to fiimax

  nr=1000
  nfii=180*2
  fiimin=40.
  fiimax=130.
  fiimin=0.
  fiimax=360.
  if(keyword_set(ff)) then begin
      fiimin=ff(0)
      fiimax=ff(1)
  endif
  if(keyword_set(nrfii)) then begin
      nr=nrfii(0)
      nfii=nrfii(1)
  endif
  r1=rr(0)
  r2=rr(1)
  
;this defines the region of the rectified image
  inter=[r1,r2,nr,nfii,fiimin,fiimax]
  
  fii=(dindgen(nfii+1)/nfii*(fiimax-fiimin)+fiimin)*!dpi/180.d0
  rad=(dindgen(nr+1)/nr*(r2-r1)+r1)
  mask=fltarr(nfii+1,nr+1)


;cassini_proje2005_f calculates original image cordinates xs,ys
;for each of the x,y,z locations defined via inter

  cassini_proje2005_f,x,y,z,xc_cas,yc_cas,b_cas,r_cas,xs,ys,$
    inter=inter,rot=rot_cas,usewa=usewa

;shift according to the pivot point
  xint=ximage0+(xs-xs0)
  yint=yimage0+(ys-ys0)

;get adu values in these original image locations via bilinear interpolation
  image_bilin,ima,xint,yint,adut
  mask(0:nfii,0:nr)=(adut)

  title=strjoin([strtrim(p),string(xs0),string(ys0)])
  if(not keyword_set(nonew)) then nwin,/large

;plot,
;together with vgr iss lit-side plot
  imaout=transpose(mask)
  mask=bytscl(imaout)
  if(keyword_set(iss)) then begin
      show_iss,[r1,r2],xs=nr+1,ys=10,lit=lit,unlit=unlit,/nop
      mask(0:nr,0:9)=bytscl(lit)
  endif
  tvplot2,mask,yy=fii*!radeg,xx=rad,/nonew,/noerase,title=title
  xima=rad
  yima=fii*!radeg
  


;--------------------------------------------------------------------------  
;stop


  if(keyword_set(encke)) then begin
      
      exclude=[153,205,250,265]
      
      rfit=133.59d0
      rmin_fit=133.5d0-1.
      rmax_fit=133.5d0+1.
      fiimin_fit=0.d0
      fiimax_fit=360.d0
      r1=rmin_fit
      r2=rmax_fit
      fiimin=fiimin_fit
      fiimax=fiimax_fit
      fiid=(dindgen(nfii+1)/nfii*(fiimax-fiimin+20.d0)+fiimin-10.d0)
      
      
;specified limits for fitting region
      ind_fii_use=where(fiid ge fiimin and fiid le fiimax,count)
      
;exclude planet, shadow ?
      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)
    indu=where(rad gt rmin_fit and rad lt rmax_fit)
    masku=mask(indu,*)
    
    for ix=0,count-1 do begin
        iix=ind_fii_use(ix)
        indmin=where(masku(*,iix) eq min(masku(*,iix)))
        rminloc(ix)=rad(indu(indmin(0)))
    endfor
    oplot,rminloc,fiid(ind_fii_use),psym=3,thick=3
    
    
    
    

    vast=''
    read,vast
    mask2=mask
    rscale=rad(1)-rad(0)
    
    for ix=0,count-1 do begin
        iix=ind_fii_use(ix)
        mask2(*,iix)=shift(mask2(*,iix),-(rminloc(ix)-rfit)/rscale)
    endfor
    tvplot2,mask2,yy=fii*!radeg,xx=rad,/nonew,/noerase,title=title        
    
    
endif


return
end