pro galfit_display_new,infile0,profile=profile,residual=residual,stop=stop,$
                       flux=flux,xr0=xr0,yr0=yr0,structure0=structure0,$
                       restore=restore,storefits=storefits,xlog0=xlog0,sample0=sample0,$
                       pa0=pa0,majorminor=majorminor,silent=silent,arcsec=arcsec,run=run,cut=cut,$
                       title0=title0,ps=ps,addps=addps,fixpos=fixpos
  
  if(n_params() le 0) then begin
     print,'============================================================================'
     print,' galfit_display_new,INFILE'
     print," visualize models corresponding to Chien Peng's GALFIT parameter file INFILE"
     print,'============================================================================'
     print,' method: "galfit -o3" -> makes "subcomps.fits" containing model components ' 
     print,'          store data to IDL-savefile "INFILE.subcomps.save" ' 
     print,'          use the data to make the desired plots'
     print,'----------------------------------------------------------------------------'
     print,' PLOT KEYWORDS:'
     print,'--------------'
     print,'  PROFILE=1       --> plot mag vs rsky profiles (plot every pixel)'
     print,'        =-1       --> improved plot (resample images before plotting)'
     print,'                          SAMPLE=sampling factor (bigger-less points)'
     print,'        =-2       --> improved plot + header information on screen'
     print,'                          SAMPLE=sampling factor (bigger-less points)'
     print,'                          TITLE=string -> print this title to the plot'
     print,'-------------'
     print,'  /RESIDUAL       --> display data-cube: obs, model, residual'
     print,'-------------'
     print,'  /MAJORMINOR     --> plot x (PA=-90) and y (PA=0) total vs obs profiles'
     print,'  = phi               plot major and minor axis profiles with major_PA=phi'
     print,'  = 999               plot major and minor axis profiles with major_PA=disk_PA'
     print,'-------------'
     print,'  /PA             --> plot major axis profiles for each component'
     print,'  = phi           --> plot major axis profiles with PA=phi'
     print,'  = 999           --> plot major and minor axis profiles with major_PA=disk_PA'
     print,'---------------------------------------------------------------------'
     print,' Keywords affecting profile plots:'
     print,'  /FLUX makes  flux vs r_sky profiles (instead of mag vs r_sky)'
     print,'  /XLOG -> use log r_sky rather than linear'
     print,'   XR=[xrmin,xrmax] and YR=[yrmin,yrmax] -> define explicit x,y range'
     print,'  /ARCSEC -> use arcsec rather than pix scale in plots'
     print,'---------------------------------------------------------'
     print,' Additional keywords:'
     print,'  /RESTORE                  --> use previously saved INFILE.subcomps.save'
     print,'                                FASTER THAN MAKING A NEW FITS-FILE'
     print,'   STRUCTURE=[st1,st2,...]  --> override structure-component labels from INFILE'
     print,'  /SILENT                   --> less output'
     print,'  /RUN                      --> call galfit_run_new,/new to make galfit decomposition' 
     print,'   CUT = npixels            --> cut image to +/- npixels around the center'
     print,'                                NOTE: the INGAL-file is modified' 
     print,'  /PS    --> make ps and png-files instead of displaying on screen'
     print,'             works with keywords:   profile=-2  --> INFILE_profile.png'
     print,'                                    residual    --> INFILE_residual.png'
     print,'                                    majorminor  --> INFILE_majorminor.png'
     print,'                                    pa  --> INFILE_pa.png'
     print,'  ADDPS=string  --> add this string to produced PS and PNG-filenames '
     print,'                    (e.g. INFILE_ADDPS_profile.png)'
     print,'---------------------------------------------------------'
     print,' VARIABLES SAVED to "INFILE.subcomps.save" ' 
     print,'      model,obs,comps,rad,ncomps,string_comps,frac_comps,zeropoint,$' 
     print,'      x1,x2,y1,y2,xc,yc,nx,ny,$' 
     print,'      data_image,out_image,sigma_image,psf_image,mask_image,xc_idl,yc_idl,$' 
     print,'      subcomp_header,obs_mask'
     print,' (see galfidl_manual.pdf)'
     print,'---------------------------------------------------------'
     print,' EXAMPLE:'
     print,"   galfit_display_new,'ngc1097_example.ingal',profile=-2,/residual,/restore"
     print,'---------------------------------------------------------'
     print,' heikki.salo@oulu.fi 220909/070110/040810'
     print,'============================================================================'
     return
  endif
  
  info=1
  if(keyword_set(silent)) then info=0

  

;----------------------------------------------------
;make ne iteration instead of displaying

  if(keyword_set(run)) then begin
     infile=remove_tag(infile0,tag='.ingal')
     galfit_run_new,/new,infile
     return
  endif

;----------------------------------------------------
;check that INFILE exist 

  close,22
  openr,22,infile0,err=stat
  if(stat ne 0) then begin
     print,infile0,' does not exist -> exit galfit_display'
     return
  endif

;--------------------------------------------------------
;check that INFILE is not an empty output, 
;produced by a crashed galfit iteration

  line=''
  readf,22,line
  if(strpos(line,'crashed') ne -1) then begin
     print,infile0,' is empty(galfit crashed) -> exit galfit_display'
     return
  endif
  close,22

;--------------------------------------------------------
;set default type/range for profile plots

  if(not keyword_set(flux)) then begin
     ytitle='magnitude/arcsec^2'
     yr=[30,15]
     if(keyword_set(xlog0)) then yr=[30,10]
  endif

  if(keyword_set(flux)) then begin
     ytitle='Flux MJy/s'
     yr=[0.01,1000]
  endif

  if(keyword_set(yr0)) then yr=yr0

;xrange defined after infile is read  
  xlog=0
  if(keyword_set(xlog0)) then xlog=xlog0


;------------------------------------------------------------
;related to optional cutting of the fit region
  cut_done=0
try_again:

;------------------------------------------------------------
; 1) read the INFILE file
;------------------------------------------------------------
; this can be ngcXXXX.ingal or outgal
; IF it contains added STRUCTURE keywords
; --> this information is displayed in the plot 
;     otherwise, comp=1, comp=2 etc are used
;------------------------------------------------------------

  input_file=infile0
  ide=infile0
  close,1
  openr,1,input_file
  line=''

;the following data will be read from infile

  data_image='unknown'
  out_image='unknown' 
  sigma_image='none'
  psf_image='none'
  mask_image='none'
  xc=-1.
  yc=-1.
  pix=0.75
  zeropoint=20.472-5*alog10(pix)
  string_comps=' '
  pix=1.
  
  mui_pos=-1                     
  mui_disk_pa=-1
  pa_disk=-90

  while not eof(1) do begin
     readf,1,line

;searching STRUCTURE lines = optional comment lines in ingal-file
;identifying the structure component

     if(strpos(line,'STRUCTURE') ne -1 and strpos(line,'STRUCTURE') le 5) $
     then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        string_comps=[string_comps,apu(2)]
                                ; print,line
     endif


;searching standard galfit input lines
if(strpos(line,'#') ne 0) then begin    ;unless a comment file

     if(strpos(line,'A)') ne -1 and strpos(line,'A)') le 2)then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        data_image=apu(1)
     endif
     if(strpos(line,'B)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        out_image=apu(1)
     endif
     if(strpos(line,'C)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        sigma_image=apu(1)
     endif
     if(strpos(line,'D)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        psf_image=apu(1)
     endif
     if(strpos(line,'F)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        mask_image=apu(1)
     endif
     if(strpos(line,'H)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        xmin=float(apu(1))
        xmax=float(apu(2))
        ymin=float(apu(3))
        ymax=float(apu(4))
        region_line=line
     endif
     if(strpos(line,'J)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        zeropoint=float(apu(1))
     endif
     if(strpos(line,'K)') ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        pix=float(apu(1))
     endif
     
;searching center position of the FIRST component,
;used in all profile plots

     npos=strpos(line,'1)')
     npos2=strpos(line,'osition')
     if(mui_pos eq -1 and npos ne -1 and npos le 2 and npos2 ne -1) then begin
        apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
        xc=float(apu(1))
        yc=float(apu(2))
        mui_pos=1
     endif

;searching disk_pa: first 'Position angle' line after 'expdisk'

     if(mui_disk_pa eq 1) then begin
        npos3=strpos(line,'osition angle')
        if(npos3 ne -1) then begin
           apu=str_sep(line,' ') & apu=apu(where(apu ne ''))
           pa_disk=float(apu(1))
           mui_disk_pa=2
        endif
     endif


;     print,line
;     print,mui_disk_pa
;     print,pa_disk
;     vast=''
;     read,vast


     if(mui_disk_pa eq -1) then begin
        if((strpos(line,'expdisk') ne -1 or strpos(line,'egdedisk') ne -1) $
           and strpos(line,'sersic') eq -1) then mui_disk_pa=1
     endif

  endif  ; not a comment

  endwhile
  close,1

;cutting
  If(keyword_set(cut) and cut_done eq 0) then begin
     ima=readfits(data_image,header,/sil)
     nx=sxpar(header,'naxis1')  
     ny=sxpar(header,'naxis2')
     xcut=cut
     ycut=cut
     xmin_new=xc-xcut
     xmax_new=xc+xcut
     ymin_new=yc-ycut
     ymax_new=yc+ycut
     i4='(i4)'
     
     sxmin=string(nint(xmin_new)>1,i4)
     symin=string(nint(ymin_new)>1,i4)
     sxmax=string(nint(xmax_new)<nx,i4)
     symax=string(nint(ymax_new)<ny,i4)
     sregion=' '+sxmin+'   '+sxmax+'  '+symin+'  '+symax
     sregion='H)  '+sregion+' # Image region to fit (xmin xmax ymin ymax)'
     job_edit2,input_file,'copy',region_line,sregion,/sil
     cut_done=1     
     goto,try_again
  endif


;default xrange
  xr_min=1
  xr_max=max([xc-xmin,yc-ymin,xmax-xc,ymax-yc])*1.4
  xr=[xr_min,xr_max]
  
;use pixels or arcsec?
  pix_use=1.
  xunit='[pix]'
  if(keyword_set(arcsec)) then begin
     pix_use=pix
     xunit='[arcsec]'
  endif
  xr=xr*pix_use

;explicitely set
  if(keyword_set(xr0)) then xr=xr0





  if(info ge 1) then begin
     print,'----------------------------------------'
     print,'Read GALFIT input-file  = ',input_file
     print,'----------------------------------------'
     print,'data_image  = ',data_image
     print,'out_image   = ',out_image
     print,'sigma_image = ',sigma_image
     print,'psf_image   = ',psf_image
     print,'mask_image  = ',mask_image
     print,'----------------------------------------'
     print,'XC,YC:',xc,yc
     print,'REGION',xmin,xmax,ymin,ymax
     print,'ZP:',zeropoint
     print,'pix:',pix
     print,'PA_DISK:', pa_disk
  endif
  
  if(keyword_set(structure0)) then begin
     print,'STRUCTURE names added via call keyword'
     string_comps=structure0
  endif

  ncomp1=n_elements(string_comps)-1
  if(ncomp1 ge 1) then begin
     string_comps=string_comps(1:ncomp1)
  endif
  


  if(keyword_set(xr0)) then xr=xr0

;-------------------------------------------------------------
; 2) Making subcomps.fits file
;-------------------------------------------------------------

;--------------------------------------------------
; Need the subcomponent fits-file ?
; if either /profile, /residual, /pa /major is set

  need_subcomp=0

  if($
     keyword_set(profile) or $ 
     keyword_set(residual) or $
     keyword_set(pa0) or $
     keyword_set(majorminor)) then $
        need_subcomp=1
  
; check if /restore is set : in this case subcomps are read from the
; previously created IDL-savefile (this is much faster) 

  if(need_subcomp eq 1 and keyword_set(restore)) then begin
     file=ide+'_subcomps.save'
     close,22 & openr,22,file,err=stat
     if(stat eq 0) then begin
        if(info ge 1) then  print,'restoring SUBCOMPS'
        restore,file
        need_subcomp=0
     endif
  endif      


;--------------------------------------------------
; Call galfit -o3 to make the subcomponent fits-file

  if(need_subcomp eq 1) then begin
     if(info ge 1) then print,'CREATE SUBCOMPS'

;subcomponents to "subcomps.fits"

     cmd='galfit -o3 '+input_file
     spawn,cmd,res     
     fits_info,'subcomps.fits',/sil,n_ext=n_ext
     
;galfit stores masked obs image

     obs=readfits('subcomps.fits',ext=1,header,/sil)
     nx=sxpar(header,'naxis1')  
     ny=sxpar(header,'naxis2')
     
     ncomps=n_ext-1
     comps=fltarr(nx,ny,ncomps)
     frac_comps=fltarr(ncomps)
     
     if(n_elements(string_comps) lt ncomps) then begin
        string_comps=replicate(' ',ncomps)
        for i=0,ncomps-1 do begin
           string_comps(i)='Comp '+string(i+1,'(i1)')
        endfor
     endif
     
;read the components

     model=obs*0.
     for i=0,ncomps-1 do begin
        comps(*,*,i)=readfits('subcomps.fits',ext=i+2,header,/sil)
        model=model+comps(*,*,i)
        frac_comps(i)=total(comps(*,*,i))
     endfor
     frac_comps=frac_comps/total(model)    
     for i=0,ncomps-1 do begin
        print,string_comps(i),' = ',frac_comps(i)
     endfor   
     
     subcomp_header=header
     
;image section

     fitsect=sxpar(header,'fitsect')
     apu=strsplit(fitsect,'[:,]',/ext)
     x1=apu(0)
     x2=apu(1)
     y1=apu(2)
     y2=apu(3)
     
;read original unmasked obs-image    

     obs=readfits(data_image,/sil)
     obs=obs(x1-1:x2-1,y1-1:y2-1)
     

;masked image

     obs_mask=obs
     if(mask_image ne 'none' ) then begin
        mask=readfits(mask_image,/sil)
        mask=mask(x1-1:x2-1,y1-1:y2-1)
        ind=where(mask ne 0,count)
        if(count ge 1) then obs_mask(ind)=0
     endif


;create array of sky distances
     
     xc_idl=xc-x1
     yc_idl=yc-y1
     dist_circle,rad,[nx,ny],xc_idl,yc_idl
     
     if(keyword_set(stop)) then stop
     
     save,file=ide+'_subcomps.save',$
          model,obs,comps,rad,ncomps,string_comps,frac_comps,zeropoint,$
          x1,x2,y1,y2,xc,yc,nx,ny,$
          data_image,out_image,sigma_image,psf_image,mask_image,xc_idl,yc_idl,$
          subcomp_header,obs_mask

;save to INFILE_subcomps.fits?

     if(keyword_set(storefits)) then begin
        cmd=' cp -rp subcomps.fits '+ide+'_subcomps.fits'
        spawn,cmd
     endif
     
  endif                         ; need_subcomp eq 1
;-------------------------------------------------------------------


;-------------------------------------------------------------
; 2) Making profile plots
;-------------------------------------------------------------

  if(keyword_set(profile)) then begin

;remove masked points from observations (since these are not fitted)
     
     obs=obs_mask

;----------------------------------------------------------
;1 original plot: plot every pixel
; can be very slow - good for small images 
;----------------------------------------------------------
     if(profile eq 1) then begin
        !p.multi=0
        nwin_hs
        col1=1
        
        if(!d.name eq 'PS') then !p.charsize=0.7
        if(!d.name eq 'PS') then col1=0


        ylog=1
        if(not keyword_set(flux)) then begin
           model=-2.5*alog10(model/pix^2>1d-16)+zeropoint
           obs=-2.5*alog10(obs/pix^2>1d-16)+zeropoint
           comps=-2.5*alog10(comps/pix^2>1d-16)+zeropoint
           ylog=0
        endif

        
        tek_color
        colors=lindgen(ncomps)+2
        
        plot,rad*pix_use,obs,ylog=ylog,psym=3,yr=yr,col=col1,xr=xr,title=ide,xlog=xlog,/nod,xtitle='rsky '+xunit
        oplot,rad*pix_use,obs,psym=3,col=120
        oplot,rad*pix_use,model,psym=3,col=col1
        for i=0,ncomps-1 do begin
           oplot,rad*pix_use,comps(*,*,i),psym=3,col=colors(i)
        endfor

        f3='(f6.3)'
        sca=1.
        if(!d.name eq 'PS') then sca=0.6
        label_data,0.45,0.92,string_comps+' '+string(frac_comps,f3),$
                   col=colors,lines=colors*0,y_i=0.05,size=1.1*sca,len=0.1
     endif




;----------------------------------------------------------
;-1 ApJL style 
;----------------------------------------------------------

     if(profile eq -1) then begin
        
        sample=1.
        if(keyword_set(sample0)) then sample=sample0
        
       print,'max(rad)',max(rad)
        if(max(rad) le 150) then sample=1


        nwin_hs
        col1=1   
        if(!d.name eq 'PS') then !p.charsize=0.7
        if(!d.name eq 'PS') then col1=0
        loadct,0,/sil
        tek_color
        colors=lindgen(ncomps)+2

        ylog=1
        if(not keyword_set(flux)) then ylog=0
        
        plot,rad*pix_use,obs,ylog=ylog,psym=3,yr=yr,col=col1,xr=xr,$
             title=ide,xlog=xlog,/nod,$
             ys=1,xtitle='Distance along sky '+xunit,ytitle=ytitle


        if(ylog eq 1) then begin   
           yp0=10^!y.crange(0)
           yp1=10^!y.crange(1)
        endif
        if(ylog eq 0) then begin
           yp0=!y.crange(0)
           yp1=!y.crange(1)
        endif
        if(xlog eq 1) then begin   
           xp0=10^!x.crange(0)
           xp1=10^!x.crange(1)
        endif
        if(xlog eq 0) then begin
           xp0=!x.crange(0)
           xp1=!x.crange(1)
        endif

        polyfill,[xp0,xp1,xp1,xp0,xp0],[yp0,yp0,yp1,yp1,yp0],col=200
        
        radi=shrink(rad,f=sample*2)  
        obsi=shrink(obs,f=sample*2)
        modeli=shrink(model,f=sample*2)


        if(not keyword_set(flux)) then begin
           obsi=-2.5*alog10(obsi/pix^2>1d-16)+zeropoint
           modeli=-2.5*alog10(modeli/pix^2>1d-16)+zeropoint
        endif

        oplot,radi*pix_use,obsi,psym=3,col=40
        oplot,radi*pix_use,modeli,psym=3,col=1


        for i=0,ncomps-1 do begin

;bulge and disk
           if(i eq 0 or i eq 1) then begin
              radi=shrink(rad,f=4*sample)       
              compsi=shrink(comps(*,*,i),f=4*sample)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              oplot,radi*pix_use,compsi,psym=3,col=colors(i)
           endif

           if(i eq 2 or string_comps(i) eq 'BAR' ) then begin
              radi=shrink(rad,f=sample)
              compsi=shrink(comps(*,*,i),f=sample)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              oplot,radi*pix_use,compsi,psym=3,col=colors(i)
           endif
           
           if(i gt 2 or string_comps(i) eq 'NUCLEUS') then begin
              radi=shrink(rad,f=0.5)
              compsi=shrink(comps(*,*,i),f=0.5)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              tabulate_f,radi,compsi,1,1000,.1,xt,mini=mini,maxi=maxi,yt
              mini=smooth(mini,5)
              maxi=smooth(maxi,5)

              npo=n_elements(mini)
              xt=xt(3:npo-5)
              mini=mini(3:npo-5)
              maxi=maxi(3:npo-5)

              polyfill,[xt,reverse(xt)]*pix_use,[maxi<yr(0),reverse(mini<yr(0))],col=colors(i)
           endif
           
        endfor

        
        f3='(f6.3)'
        sca=1.
        if(!d.name eq 'PS') then sca=0.6

        label_data,0.45,0.92,string_comps+' '+string(frac_comps,f3),$
                   col=colors,lines=colors*0,y_i=0.05,size=1.1*sca,len=0.1
        
        label_data,0.1,0.92,['OBS','MODEL'],$
                   col=[140,1],lines=[0,0],y_i=0.075,size=1.1*sca,len=0.1
     endif
     



;----------------------------------------------------------
; ApJL style + info
;----------------------------------------------------------

     if(profile eq -2) then begin
        
        psfile=infile0+'_profile'
        if(keyword_set(addps)) then psfile=infile0+addps+'_profile'
        if(keyword_set(ps)) then psdirect,psfile,-2,/color

        sample=1.
        if(keyword_set(sample0)) then sample=sample0

        print,'max(rad)',max(rad)
        if(max(rad) le 150) then sample=1


        if(not keyword_set(fixpos)) then nwin_hs
        if(    keyword_set(fixpos)) then nwin_hs,xpos=20,ypos=50

        col1=1   
        if(!d.name eq 'PS') then !p.charsize=0.7
        if(!d.name eq 'PS') then col1=0
        loadct,0,/sil
        tek_color
        colors=lindgen(ncomps)+2

        ylog=1
        if(not keyword_set(flux)) then ylog=0
        
        plot,rad*pix_use,obs,ylog=ylog,psym=3,yr=yr,col=col1,xr=xr,$
             title=ide,xlog=xlog,/nod,$
             ys=1,xtitle='Distance along sky '+xunit,ytitle=ytitle,position=[0.1,0.1,0.75,0.75]

        if(ylog eq 1) then begin   
           yp0=10^!y.crange(0)
           yp1=10^!y.crange(1)
        endif
        if(ylog eq 0) then begin
           yp0=!y.crange(0)
           yp1=!y.crange(1)
        endif
        if(xlog eq 1) then begin   
           xp0=10^!x.crange(0)
           xp1=10^!x.crange(1)
        endif
        if(xlog eq 0) then begin
           xp0=!x.crange(0)
           xp1=!x.crange(1)
        endif

        polyfill,[xp0,xp1,xp1,xp0,xp0],[yp0,yp0,yp1,yp1,yp0],col=200
        
        radi=shrink(rad,f=sample*2)  
        obsi=shrink(obs,f=sample*2)
        modeli=shrink(model,f=sample*2)


        if(not keyword_set(flux)) then begin
           obsi=-2.5*alog10(obsi/pix^2>1d-16)+zeropoint
           modeli=-2.5*alog10(modeli/pix^2>1d-16)+zeropoint
        endif

        oplot,radi*pix_use,obsi,psym=3,col=80
        oplot,radi*pix_use,modeli,psym=3,col=1


;to make smaller plots:
; sample bulge by 4*sample
;        disk  by 4*sample
;        bar   by   sample
; other component = fill with polygons



        for i=0,ncomps-1 do begin


           compsi=-2.5*alog10(comps(*,*,i)/pix^2>1d-16)+zeropoint
           if(i eq 0) then print,'npix',n_elements(compsi)
           junk=where(compsi lt 30,junkcount)
           print,string_comps(i),junkcount

           if(i eq 0) then begin



              radi=rad
              compsi=comps(*,*,i)

              f=junkcount/2000.
              if(f gt 1.5) then begin
                 f=nint(f)<4

              radi=shrink(rad,f=4*sample)       
              compsi=shrink(comps(*,*,i),f=4*sample)
              endif

              compsi_mag=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              junki=where(compsi_mag lt 30,junkcount)

              if(not keyword_set(flux)) then compsi=compsi_mag
              if(junkcount ge 2) then $
                 oplot,radi(junki)*pix_use,compsi(junki),psym=3,col=colors(i)
           endif


           if(i eq 1) then begin
              radi=shrink(rad,f=4*sample)       
              compsi=shrink(comps(*,*,i),f=4*sample)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              oplot,radi*pix_use,compsi,psym=3,col=colors(i)
           endif

           if(i eq 2 or string_comps(i) eq 'BAR' ) then begin
              radi=shrink(rad,f=sample)
              compsi=shrink(comps(*,*,i),f=sample)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              oplot,radi*pix_use,compsi,psym=3,col=colors(i)
           endif
           
           if(i gt 2 or string_comps(i) eq 'NUCLEUS') then begin
              radi=shrink(rad,f=0.5)
              compsi=shrink(comps(*,*,i),f=0.5)
              if(not keyword_set(flux)) then $
                 compsi=-2.5*alog10(compsi/pix^2>1d-16)+zeropoint
              tabulate_f,radi,compsi,1,1000,.1,xt,mini=mini,maxi=maxi,yt
              mini=smooth(mini,5)
              maxi=smooth(maxi,5)

              npo=n_elements(mini)
              xt=xt(3:npo-5)
              mini=mini(3:npo-5)
              maxi=maxi(3:npo-5)
              polyfill,[xt,reverse(xt)]*pix_use,[maxi<yr(0),reverse(mini<yr(0))],col=colors(i)
           endif
           
        endfor

        f3='(f6.3)'
        sca=1.
        if(!d.name eq 'PS') then sca=0.6

        label_data,0.45,0.92,string_comps+' '+string(frac_comps,f3),$
                   col=colors,lines=colors*0,y_i=0.05,size=1.1*sca,len=0.1,thick=colors*0+3
        
        label_data,0.1,0.92,['OBS','MODEL'],$
                   col=[140,1],lines=[0,0],y_i=0.075,size=1.1*sca,len=0.1,thick=[3,3]


;add info from subcomp_header
        plot,lindgen(2),pos=[0.,0.0,0.97,0.97],/nod,/noerase,XS=15,YS=15
        dy=0.025
        ypos=1.0

        for i=7,12 do begin
           line=Subcomp_HEADER(i)

           if(strpos(line,'COMMENT') ne -1) then line='-----------------------'
           posi=strpos(line,'/')
           if(n_elements(posi) ne 0) then begin
              if(n_elements(posi) eq 1) then begin
                 line=strmid(line,0,posi-1)
              endif
              if(n_elements(posi) eq 2) then begin
                 line=strmid(line,0,posi(1)-1)
              endif
           endif

           ypos=ypos-dy
           xyouts,0.05,ypos,line,/data
        endfor

        ypos=1.0
        for i=13,15 do begin
           line=Subcomp_HEADER(i)

           if(strpos(line,'COMMENT') ne -1) then line='-------------------------'
           posi=strpos(line,'/')
           if(n_elements(posi) ne 0) then begin
              if(n_elements(posi) eq 1) then begin
                 line=strmid(line,0,posi-1)
              endif
              if(n_elements(posi) eq 2) then begin
                 line=strmid(line,0,posi(1)-1)
              endif
           endif           
           ypos=ypos-dy
           xyouts,0.45,ypos,line,/data
        endfor

        ypos=ypos-dy
        xyouts,0.45,ypos,!stime,/data
        ypos=ypos-dy
        xyouts,0.45,ypos,mywhoami()+'@'+hostname(),/data

        ypos=ypos-dy
        xyouts,0.45,ypos,/data,'code version: heikki.salo@oulu.fi 290710',chars=sca*.8



        plot,lindgen(2),pos=[0.78,0.05,0.97,0.97],/nod,/noerase,XS=15,YS=15
        nlines=n_elements(subcomp_HEADER)
        dy=0.93/(nlines-2*ncomps-25)<0.04
        ypos=1.0

        if(keyword_set(title0)) then begin 
           xyouts,0.01,ypos,title0
           ypos=ypos-dy
        endif

        mui=0
        okmui=0
        for i=17,nlines-1 do begin
           line=Subcomp_HEADER(i)
          line=strjoin(strsplit(line,"'",/ext))
 ;          if(strpos(line,'Degrees of') ne -1) then mui=1
          if(strpos(line,'CHI') ne -1) then mui=1

           ok=1
           if(strpos(line,'_XC') ne -1) then ok=0
           if(strpos(line,'_YC') ne -1) then ok=0

           if(strpos(line,'COMMENT') ne -1) then line='-------------------------'
           posi=strpos(line,'/')
           if(n_elements(posi) ne 0) then begin
              if(n_elements(posi) eq 1) then begin
                 line=strmid(line,0,posi-1)
              endif
              if(n_elements(posi) eq 2) then begin
                 line=strmid(line,0,posi(1)-1)
              endif
           endif

           step=1.
           if(strtrim(line,2) eq '') then step=.5

           if(mui eq 0 and ok eq 1) then begin
              ypos=ypos-dy*step
              xyouts,0.01,ypos,line,/data
              print,line
           endif
        endfor


;calculate chisq here

     res=obs-model
     obs_mask=obs
     res_mask=res
     npix=n_elements(res)
     if(mask_image ne 'none' ) then begin
        mask=readfits(mask_image,/sil)
        mask=mask(x1-1:x2-1,y1-1:y2-1)
        ind=where(mask ne 0,count)
        npix=npix-count
        if(count ge 1) then obs_mask(ind)=0
        res_mask=res
        if(count ge 1) then res_mask(ind)=0
     endif

        sigma=readfits(sigma_image,/sil)
        sigma=sigma(x1-1:x2-1,y1-1:y2-1)
        chi2_title=total(res_mask^2/sigma^2)/npix
        xyouts,0.01,ypos-dy,/data,'Chi2/nu='+string(chi2_title)

;        endif

        if(keyword_set(ps)) then psdirect,psfile,-2,/color,/stop
     endif                         


  endif                         ;end /profile 


;-------------------------------------------------------------
; 4) Make residual plots ?
;-------------------------------------------------------------

  if(keyword_set(residual)) then begin

     psfile=infile0+'_residual'
     if(keyword_set(addps)) then psfile=infile0+addps+'_residual'
     if(keyword_set(ps)) then psdirect,psfile,-2,/color


     if(info ge 1) then       print,'restoring subcomps.save'
     restore,ide+'_subcomps.save'


     res=obs-model
     obs_mask=obs
     res_mask=res

     npix=n_elements(res)

     if(mask_image ne 'none' ) then begin
        mask=readfits(mask_image,/sil)
        mask=mask(x1-1:x2-1,y1-1:y2-1)
        ind=where(mask ne 0,count)
        npix=npix-count
        if(count ge 1) then obs_mask(ind)=0
        res_mask=res
        if(count ge 1) then res_mask(ind)=0
     endif

        sigma=readfits(sigma_image,/sil)
        sigma=sigma(x1-1:x2-1,y1-1:y2-1)

        chi2_title=total(res_mask^2/sigma^2)/npix
        print,'Chi2/nu=', chi2_title


;stop

     loadct,0,/sil
     !p.multi=[0,2,2]
     xr1=[1,1000]
     yr1=[0.01,1000]

     scale=1.*(float(y2)-float(y1))/(float(x2)-float(x1))


        if(not keyword_set(fixpos)) then nwin_hs,xs=800,ys=800*scale
        if(    keyword_set(fixpos)) then nwin_hs,xpos=720,ypos=50,xs=800,ys=800*scale


     nx=n_elements(obs(*,0))
     ny=n_elements(obs(0,*))
     xx=lindgen(nx)-xc+x1
     yy=lindgen(ny)-yc+y1

     xx=xx*pix_use
     yy=yy*pix_use

     tvplot2,alog(obs>yr1(0)<yr1(1)),xx=xx,yy=yy,/asp,title=ide ,/nonew,xtitle=xunit,ytitle=xunit
     tvplot2,alog(obs_mask>yr1(0)<yr1(1)),xx=xx,yy=yy,/asp,title=ide+' masked',xtitle=xunit,ytitle=xunit
     tvplot2,alog(model>yr1(0)<yr1(1)),xx=xx,yy=yy,/asp,/nonew,title='model',xtitle=xunit,ytitle=xunit
     tvplot2,res_mask>(-10*yr1(0))<yr1(0)*10,xx=xx,yy=yy,/asp,/nonew,title='residual'+' '+string(chi2_title),xtitle=xunit,ytitle=xunit
     !p.multi=0
     if(keyword_set(ps)) then psdirect,psfile,-2,/color,/stop

  endif


;-------------------------------------------------------------
; 5) Make major/minor plots ?
;-------------------------------------------------------------


  if(keyword_set(majorminor)) then begin

     psfile=infile0+'_majorminor'
     if(keyword_set(addps)) then psfile=infile0+addps+'_majorminor'
     if(keyword_set(ps)) then psdirect,psfile,-2,/color


     restore,ide+'_subcomps.save'
     obs=obs_mask

;major axis phi
; /major         --> use PA=-90
;  major=-999    --> use PA=disk PA

     pa=majorminor
     if(pa eq 1) then pa=-90
     if(pa eq 999) then pa=pa_disk
     phi=(pa+90.)/!radeg        ;angle from x-axis in radians

     nrad=min([xc-float(x1),float(x2)-xc,yc-float(y1),float(y2)-yc])-1
     radi=findgen(nrad*2+1)-nrad
     
     xx=xc_idl+radi*cos(phi)
     yy=yc_idl+radi*sin(phi)
     image_bilin,obs,xx,yy,obs_major
     image_bilin,model,xx,yy,model_major
     
     xx=xc_idl+radi*cos(phi+!pi/2)
     yy=yc_idl+radi*sin(phi+!pi/2)
     image_bilin,obs,xx,yy,obs_minor
     image_bilin,model,xx,yy,model_minor

     tek_color
     nwin_hs
     yr=[30,10]
     if(keyword_set(yr0)) then yr=yr0

     plot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,/nod,$
          yr=yr,title=ide,xtitle='radius '+xunit,ytitle='mag'
     oplot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,col=2
     oplot,radi*pix_use,-2.5*alog10(obs_minor/pix^2>1d-16)+zeropoint,col=3       
     oplot,radi*pix_use,-2.5*alog10(model_major/pix^2>1d-16)+zeropoint,col=8,lines=2
     oplot,radi*pix_use,-2.5*alog10(model_minor/pix^2>1d-16)+zeropoint,col=9,lines=2
     
     label_data,0.1,0.9,['major pa='+string(pa,'(f5.1)'),'minor pa='+string(pa+90,'(f5.1)')],col=[2,3]

;small insert plot

     obs_major=(obs_major+reverse(obs_major))/2.
     obs_minor=(obs_minor+reverse(obs_minor))/2.
     plot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,/nod,yr=yr,/xlog,$
          xr=xr,title='',/noe,pos=[0.63,0.7,0.93,0.93]
     oplot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,col=2
     oplot,radi*pix_use,-2.5*alog10(obs_minor/pix^2>1d-16)+zeropoint,col=3       
     oplot,radi*pix_use,-2.5*alog10(model_major/pix^2>1d-16)+zeropoint,col=8,lines=2
     oplot,radi*pix_use,-2.5*alog10(model_minor/pix^2>1d-16)+zeropoint,col=9,lines=2
     
     label_data,0.1,0.9,['major','minor'],col=[2,3],lines=[0,0]

     if(keyword_set(ps)) then psdirect,psfile,-2,/color,/stop
     
  endif


;-------------------------------------------------------------
; 5) Make phi=PA plots for each component?
;-------------------------------------------------------------

  
  if(keyword_set(pa0)) then begin

     psfile=infile0+'_pa'
     if(keyword_set(addps)) then psfile=infile0+addps+'_pa'
     if(keyword_set(ps)) then psdirect,psfile,-2,/color

     restore,ide+'_subcomps.save'
     obs=obs_mask

;major axis phi
; /pa         --> use PA=-90
;  pa=999     --> use PA=disk PA

     pa=pa0
     if(pa eq 1) then pa=-90
     if(pa eq 999) then pa=pa_disk
     phi=(pa+90.)/!radeg        ;angle from x-axis in radians

     nrad=min([xc-float(x1),float(x2)-xc,yc-float(y1),float(y2)-yc])-1
     radi=findgen(nrad*2+1)-nrad
     nrad=n_elements(radi)


     col1=1        
     if(!d.name eq 'PS') then !p.charsize=0.7
     if(!d.name eq 'PS') then col1=0
     

     obs=obs/pix^2
     model=model/pix^2
     comps=comps/pix^2

     xx=xc_idl+radi*cos(phi)
     yy=yc_idl+radi*sin(phi)
     image_bilin,obs,xx,yy,obs_major
     image_bilin,model,xx,yy,model_major
     
     xx=xc_idl+radi*cos(phi+!pi/2)
     yy=yc_idl+radi*sin(phi+!pi/2)
     image_bilin,obs,xx,yy,obs_minor
     image_bilin,model,xx,yy,model_minor
     
     apu=findgen(nrad)
     comps_major=findgen(nrad,ncomps)
     comps_minor=findgen(nrad,ncomps)
     for i=0,ncomps-1 do begin
        xx=xc_idl+radi*cos(phi)
        yy=yc_idl+radi*sin(phi)
        image_bilin,comps(*,*,i),xx,yy,apu
        comps_major(*,i)=apu

        xx=xc_idl+radi*cos(phi+!pi/2)
        yy=yc_idl+radi*sin(phi+!pi/2)
        image_bilin,comps(*,*,i),xx,yy,apu
        comps_minor(*,i)=apu
     endfor

     yr=[30,10]
     if(keyword_set(yr0)) then yr=yr0
     nwin_hs
     tek_color
     colors=lindgen(ncomps)+2

     obs_major=(obs_major+reverse(obs_major))/2.
     obs_minor=(obs_minor+reverse(obs_minor))/2.
     
     plot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,/nod,$
          yr=yr,title=ide,xtitle='radius '+xunit,$
          ytitle='mag',xr=xr,/xlog,ns=3
     
     oplot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,thick=3,lines=2,ns=3
     
     oplot,radi*pix_use,-2.5*alog10(model_major/pix^2>1d-16)+zeropoint
     label_data,0.1,0.94,['OBS','MODEL'],lines=[2,0],col=[col1,col1],thick=[4,2]


     for i=0,ncomps-1 do begin      
        oplot,radi*pix_use,-2.5*alog10(comps_major(*,i)/pix^2>1d-16)+zeropoint,col=colors(i)
     endfor

     plot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,/nod,yr=yr,xr=xr,$
          title='',/noe,pos=[0.63,0.63,0.93,0.9],xs=1
     oplot,radi*pix_use,-2.5*alog10(obs_major/pix^2>1d-16)+zeropoint,col=2,lines=2
     oplot,radi*pix_use,-2.5*alog10(obs_minor/pix^2>1d-16)+zeropoint,col=3,lines=2
     oplot,radi*pix_use,-2.5*alog10(model_major/pix^2>1d-16)+zeropoint,col=8,lines=0
     oplot,radi*pix_use,-2.5*alog10(model_minor/pix^2>1d-16)+zeropoint,col=9,lines=0
     label_data,0.1,0.9,['major axis','minor axis'],col=[2,3],lines=[0,0],y_i=0.2,shift=-0.2

     if(keyword_set(ps)) then psdirect,psfile,-2,/color,/stop
  endif

!p.charsize=1
end