; Copyright (C) 2016      David Tsiklauri <d.tsiklauri@qmul.ac.uk>
;
; This program is free software: you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation, either version 3 of the License, or
; any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details:
; <http://www.gnu.org/licenses/>.

ns=14 ; this many data points will be subracted from the fit
n=20
nx=20000

f=dblarr(nx+4)
fphd=dblarr(nx,n+1)

x=dblarr(nx+4)
xph=dblarr(nx)

t=dblarr(n+1)
fm=dblarr(n+1)

tl=dblarr(n+1-ns)
fl1=dblarr(n+1-ns)

tt=0.0d0

OPENR, 1, '00',  /F77_UNFORMATTED
READU, 1,f,x,tt
CLOSE, 1

for i=0,nx-1 do begin
fphd(i,0)=f(i+2)
xph(i)=x(i+2)
endfor

window,0
plot,xph,fphd(*,0),xrange=[min(xph)/10.0,max(xph)/10.0],xstyle=1,yrange=[-1,1],ystyle=1

for h=0,n do begin

if (h LE 9) then begin

OPENR, 1, '0'+strtrim(string(h),2),  /F77_UNFORMATTED
READU, 1,f,x,tt
CLOSE, 1
endif else begin 
OPENR, 1, ''+strtrim(string(h),2),  /F77_UNFORMATTED
READU, 1,f,x,tt
CLOSE, 1
endelse

for i=0,nx-1 do begin
fphd(i,h)=f(i+2)
xph(i)=x(i+2)
endfor

t(h)=tt
fm(h)=max(fphd(*,h))

print,h,t(h),max(fphd(*,h))

oplot,xph,fphd(*,h),line=1

endfor

for i=0,n-ns do tl(i)=t(i+ns)
for i=0,n-ns do fl1(i)=fm(i+ns)

r=poly_fit(alog(tl),alog(fl1),1,yfit,SIGMA=sigma)
print,r

window,1
plot,t,fm,/xlog,/ylog,xrange=[1,20],xstyle=1,yrange=[0.1,1],ystyle=1
oplot,t,exp(r[0])*t^r[1], psym=5
save,r,fphd,fm,xph,t,filename='disp_dif.sav'

end