pro density_find
common betafit, a

;d.sur_bri = cnt/pix^2
;exposure in d.exposure
; d.r in pixel

d = mrdfits('sur_bri_con.fits',1)
r = 0.5*(d.r[0]+d.r[1])
x = where(r le 1000)

;Here make sure r in cm, sur_bri in cnts/s/arcmin^2
arcsec2cm = (1.0/3600.0)*(!pi/180.)*1.11d27*0.89

rnew = r(x)*0.492*arcsec2cm

sb = d(x).sur_bri/d(x).exposure * (60.0*60.0/0.492/0.492)
sberr = d(x).sur_bri_err/d(x).exposure * (60.0*60.0/0.492/0.492)

n = n_elements(rnew)
!p.multi=[0,1,2]


;Fit this data with a beta model
weights=1.0/sberr/sberr

a = [0.1, 200.0*0.492*arcsec2cm, 0.6]
res = curvefit(rnew, sb, weights, a, function_name='beta', itmax=1000)

print, 'Surface brightness fit result:'
print, a


plot, rnew, sb, /xlog, /ylog, psym=10
oplot, rnew, res
for i=0,n-1 do oplot, [rnew(i),rnew(i)], sb(i)+[-sberr(i),sberr(i)]


;Now find counts from surface brightness
;First need to resort such that everything goes from outside in.
sb_n = reverse(sb)
sberr_n = reverse(sberr)

rin = reverse(d(x).r[0]*0.492*arcsec2cm) & rout = reverse(d(x).r[1]*0.492*arcsec2cm)

rout = [rout,0]

area = !pi*(rout^2 - rin^2)
dist = 1.11d27*0.89
b = (!pi*dist/10800)^2

c = dblarr(n)

;create volume array
V = dblarr(n+1,n+1)
for j=0, n do begin
    for i=0, n do begin
        if (j gt i) then $
          v(i,j) = 4*!pi/3 *(rout(i)^2 - rout(j)^2)^(3./2.) $
          else v(i,j) = 0
    endfor
endfor

c(0) = sb_n(0)*area(0)/b/v(0,1)

v_Sphere = 4.*!pi*rout^3/3.

for m=1,n-1 do begin
    temp = 0
    for i=0,m-1 do $
      temp = temp + c(i)*((v(i,m+1) - v(i+1,m+1)) - (v(i,m) - v(i+1,m)))
    
    c(m) = (sb_n(m)*area(m)/b - temp)/v(m,m+1)
endfor

c_total = c*v_Sphere

;need conversion from counts/s ==> apec normalization
;then from apec normalization ==> electron density

;to get conversion from counts/s ==> apec normalization
ctspers = [1.830,0.1832,0.01854,0.00185]
ctspers_err = [1.4e-3,4e-4,1.4e-4,5e-5]
apec_norm = [1e-2,1e-3,1e-4,1e-5]

conv = total(ctspers*apec_norm/ctspers_err^2)/total(ctspers^2/ctspers_err^2)
a_norms = c_total*conv
n_e = sqrt(a_norms*4*!pi*dist^2*(1+0.0943)^2/(1d-14*v_Sphere))

plot, reverse(rnew), n_e, psym=10, /xlog, /ylog

b = [0.005]
ne_fit = curvefit(reverse(rnew), n_e, none, b, function_name='nebeta', itmax=1000)
print, b, max(n_e)
oplot, reverse(rnew), ne_fit

stop

end


pro nebeta, X, B, F, pder
common betafit, a
  f = b[0]*(1 + (X/a[1])^2)^(-3*a[2]/2.)

IF N_PARAMS() GE 4 THEN $   
  pder = [[(1 + (X/a[1])^2)^(-3*a[2]/2.)]]
end


pro beta, X, A, F, pder
  f = a[0]*(1 + (X/a[1])^2)^(-3*a[2] + 0.5)

IF N_PARAMS() GE 4 THEN $   
  pder = [[(1 + (X/a[1])^2)^(-3*a[2] + 0.5)],[2*a[0]*X*(-3*a[2] + 0.5)*(1 + (X/a[1])^2)^(-3*a[2] - 0.5) /a[1]/a[1]],[-3*a[0]*(1 + (X/a[1])^2)^(-3*a[2] + 0.5)*alog(1 + (X/a[1])^2)]]

end