; IDL program to plot (and fit?) the azimuthally averaged radial 
; of profile of the GRB030329 host galaxy.
; Setting ngrid to an integer larger than 1 will cause the profil 
; to be measured with a subsampling of ngrid (takes longer time 
; of course).

!P.charsize=1.5
!P.MULTI = [0,1,2]

nsub = 900 
ngrid=1
pixscale = 1.0
;fit radii
minr = 30.
maxr =250.

;First read in image
im = readfits('NGC4736I_comb.fits')
mask = readfits('mask.fits')
nim = n_elements(im(0,*))

;Subtract background
im = im-960.
writefits,'skysub.fits',im

;Cut out a subimage
xc=940.
yc=1140.
subim = fltarr(nsub,nsub)
submask = fltarr(nsub,nsub)
for nx=0,nsub-1 do begin
   for ny=0,nsub-1 do begin
       subim(nx,ny) = im(xc-nsub/2+nx,yc-nsub/2+ny)
       submask(nx,ny) = mask(xc-nsub/2+nx,yc-nsub/2+ny)
   endfor
endfor

;Make NxN subgrid image
imsubgrid = fltarr(ngrid*nsub,ngrid*nsub)
masksubgrid = fltarr(ngrid*nsub,ngrid*nsub)
for nx=0,ngrid*nsub-1 do begin
   for ny=0,ngrid*nsub-1 do begin
      imsubgrid(nx,ny) = subim(nx/ngrid,ny/ngrid)
      masksubgrid(nx,ny) = submask(nx/ngrid,ny/ngrid)
   endfor
endfor

;Fit centre
fit = gauss2dfit(imsubgrid,A)
xcen = fix(A(4)+0.5)
ycen = fix(A(5)+0.5)
print,'Center:',xcen,ycen
ii = where(imsubgrid GE 2.^16.)
imsubgrid(ii) = 2.^16.

;Average profile in circles around this centre
;Inspired by http://www.dfanning.com/tips/where_to_2d.html
xx = shift( dist(n_elements(imsubgrid(0,*))),xcen,ycen)
d = findgen(ngrid*nsub/2)
profile = d*0.
scatter = d*0.
for n=0,n_elements(d)-2 do begin
  ii = where(xx GE d(n) and xx LT d(n+1) and masksubgrid NE 1)
  s = SIZE(xx)
  ncol = s(1)
  col = ii MOD ncol
  row = ii / ncol
  profile(n) = mean(imsubgrid(col,row))
  if (n_elements(imsubgrid(row,col)) GE 2) then scatter(n) = stddev(imsubgrid(row,col))/sqrt(n_elements(imsubgrid(row,col)))
endfor

;Get rid of 0 at the end of the array
d = d(0:n_elements(d)-2)
profile = profile(0:n_elements(profile)-2)
scatter = scatter(0:n_elements(scatter)-2)

;Plot the result
window,xsize=600,ysize=800
radius = d/float(ngrid)
plot_io,radius,profile,yrange=[0.1,max(profile)*1.2],$
xrange=[0.1,1.05*max(radius)],psym=3,/ystyle,/xstyle,$
ytitle='Surface Brightness',title='NGC4736 I-band',position=[0.15,0.32,0.95,0.95],$
xtickname=[' ',' ',' ',' ',' ',' ']
errplot,radius,profile-scatter,profile+scatter


;Fit the outer parts with a sersic function using mpfit
;http://cow.physics.wisc.edu/~craigm/idl/fitting.html
col = getcolor(/load)
ii = where(radius GE minr and radius LE maxr)
A = [13000.,13.,0.25]
fit = MPFITFUN('sersic',radius(ii),profile(ii),scatter(ii),A,COVAR=cov)
print
print,'FINAL RESULTS:'
oplot,radius(ii),sersic(radius(ii),fit),linestyle=2,color=col.red,thick=3
print,'Sersic Halflight radius:',fit[1]*pixscale,'   +/-',sqrt(cov(1,1))*pixscale,$
        ' arcsec'
print,'Sersic n:',1./fit[2],'   +/-',sqrt(cov(2,2))

;Plot residuals
plot,radius(ii),(profile(ii)-sersic(radius(ii),fit))/profile(ii)*100.,$
xrange=[0.3,1.05*max(radius)],/ystyle,/xstyle,$
xtitle='R (pixles)',ytitle='Residuals (%)',position=[0.15,0.10,0.95,0.32]





end