; 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
nk=40000
k=dblarr(nk)

f=dblarr(nx+4)
fph=dblarr(nx,n+1)
fan=dblarr(nx,n+1)

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

t=dblarr(n+1)
fm1=dblarr(n+1)
fm2=dblarr(n+1)

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

tt=0.0d0

for i=0,nk-1 do k(i)=-250.0+500.0*i/float(nk-1)

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

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

window,0
plot,xph,fph(*,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 j=0,nx-1 do begin
fph(j,h)=f(j+2)
xph(j)=x(j+2)
endfor

t(h)=tt
for i=0,nx-1 do begin
fan(i,h)=int_tabulated(k,exp(-k^2/4.0)*cos(k^3*t(h)+k*xph(i))/(2.0*sqrt(!pi)))
endfor

fm1(h)=max(fph(*,h))
fm2(h)=max(fan(*,h))

print,h,t(h),max(fan(*,h)),max(fph(*,h)), max(fan(*,h))-max(fph(*,h))

oplot,xph,fph(*,h),line=1
oplot,xph,fan(*,h),psym=4

endfor

for i=0,n-ns do tl(i)=t(i+ns)
for i=0,n-ns do fl1(i)=fm1(i+ns)
for i=0,n-ns do fl2(i)=fm2(i+ns)

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

r2=poly_fit(alog(tl),alog(fl2),1,yfit,SIGMA=sigma)
print,r2

window,1
plot,t,fm1,/xlog,/ylog,xrange=[1,20],xstyle=1,yrange=[0.25,1],ystyle=1
oplot,t,exp(r1[0])*t^r1[1], psym=5
oplot,t,exp(r2[0])*t^r2[1], psym=1
save,r1,r2,fan,fph,filename='disp_kdv.sav'

end