pro theory,  day, tmin, tmax, ps_or_x
; modified from theory_stand_alone.pro 4/4/95
 
;theory.pro in ~pete/de/ICRH/theory/theory.pro
;
; read wave and particle distribution data scaled from
; 83064, 82047, and 82111 data files for three masses
; use the fw velocity numbers to calculate characteristic energies
; and compare those with the asymptotic value inferred from
; the power spectral density and slope of powerspectral density
; using Crew's asymptiotic formulation.....
;	wkp 3/95
;  uses the poly_fit idl routine.
yyddd=''
title='' 
mass_name=strarr(3)
mass_name=['he','o','h']
mass=fltarr(3)
mass=[4,16,1]
;print,mass_name,mass

dum=fltarr(9)		; 9 values on a line
sec=fltarr(10)		; 10 max time intervals in file
hour=fltarr(10)		;  time in decimal hours
re=fltarr(10)		; geocentric distance in re for each time
hz=fltarr(10)		; gyrofreq at each time in hz
alpha=fltarr(10)	; alpha parameter for each time
psd=fltarr(10)		; psd in v^2/m^2/s for each time at gyrofreq
vel=fltarr(4)		; full width at 4 selected fraction of Fmax
echar=fltarr(10)	; characteristic energy derived using poly_fit
;			; from the (up to) 4 velocity values
etheory=fltarr(10)	; characterictic energy using eqn 13 from Retterer
fudge=fltarr(10) 	; echar/etheory   
NDegree=1L		; use a first order poly_fit
ff=fltarr(4)		; set fractions of Fmax for which vel are given
ff(0)=alog(0.3162278)
ff(1)=alog(0.1000000)
ff(2)=alog(0.0316228)
ff(3)=alog(0.0100000)
;print,'FF: ',ff
; factors for Retter's eqn 13.
q=1.0
qp=1.6e-19
mp=1.67e-27
eV=1.6e-19
ReR=6367e3

for Mnum=0,2 do begin
; loop throough masses, but start with just he
;
!P.multi=0

file_name=mass_name(Mnum)+string(day,format='(i5)')+'.dat'
mass_value=mass(Mnum)
print,file_name
openr,getlun,file_name,/get_lun
lun=getlun
END_OF_FILE=0
; read 9 lines the last is the plot title
for ii=0,8 do begin
 readf,lun,title
 print,title
endfor
ii=-1
; loop through the xx values for each mass
; start with just one value

repeat begin				;BEGIN REPEAT
if eof(lun) then begin
 END_OF_FILE=1
 print, 'END OF FILE'
 goto, end_of_mass
endif
ii=ii+1
readf,lun,dum
;print,dum
sec(ii)=dum(0)
hour(ii)=dum(0)/3600.
re(ii)=dum(1)
hz(ii)=dum(2)
alpha(ii)=dum(3)
psd(ii)=dum(4)
for i=0,3 do vel(i)=dum(5+i)/2.0	; half velocity is relevant!!!
;print, 'HALF vel',vel
;echar=fltarr(10)	; characteristic energy derived using poly_fit
;			; from the (up to) 4 velocity values
;ff=fltarr(10)		; set fractions of Fmax for which vel are given
; some values are 0 .... ignore them for the fit.
kk=3			; max of 4 velocity values

if (vel(0) eq 0.0 ) then kk=2
if (vel(1) eq 0.0 ) then kk=1
if (vel(2) eq 0.0 ) then begin
 print, 'Not enough velocity values ',vel
 kk=0
endif
; make polyfit
print,'sec/number of points/vel',sec(ii),kk+1,vel

E=fltarr(kk+1)			; fit 4,3,or 2 points
F=fltarr(kk+1)

if( kk gt 0 ) then begin
 k1=-1
 for k= 3-kk,3 do begin
  k1=k1+1
  e(k1)=mass(Mnum)*(vel(k)/13.84)^2
  f(k1)=ff(k)
 endfor 
endif
print,'e and f',e,f
result=poly_fit(e,f,NDegree)
;print,'RESULT=',result
if result(1) eq 0.0 then begin
	print, 'ECHAR=0 TIME=  ',file_name,sec(ii)
	echar(ii)=0.0
endif else begin
	echar(ii)=-1.0/result(1)
endelse
; calculate the theortical characteristic energy
eta=1.0
print, 'You may want to vary eta by mass '
aa=alpha(ii)
pp=psd(ii)
mm=mass_value
ss=re(ii)
se=psd(ii)
dd=se*0.25*eta*(q*qp/(mm*mp))^2
print, 'Diffusion coeff',dd
v0=(ss*ReR*dd)^0.333333
print, 'v0 (km/s) =',v0
af=(3.0*aa+5.5)^0.333333 / (3.*aa+1.0)^0.666667
print, 'af=',af
WW=af*mm*mp*v0^2
WW=WW/eV
print, 'Theoretical Characteristic Energy (eV)',WW
etheory(ii)=WW	
if echar(ii) ne 0.0 then begin
	fudge(ii)=etheory(ii)/echar(ii)
endif else begin
	fudge(ii)=0.0
endelse



print, 'END of one time interval increment to next'
endrep until (END_OF_FILE eq 1)
end_of_mass:
close,lun


print,echar
xstyle_num=0

plot,hour(0:ii),echar(0:ii),title=title,ytitle='Characteristic Energy (eV)', $
 xtitle='Time in decimal Hours',xrange=[tmin/3600.,tmax/3600.], $
position=[0.15,0.52,0.85,0.72],xstyle=1,/normal

plot,hour(0:ii),etheory(0:ii),title=title,ytitle='Theoreitical  Energy (eV)', $
xtitle='Time in decimal Hours',xrange=[tmin/3600.,tmax/3600.], $
position=[0.15,0.41,0.85,0.61],xstyle=1,/normal

plot,hour(0:ii),fudge(0:ii),title=title,ytitle='Correction factor', $
xtitle='Time in decimal Hours',xrange=[tmin/3600.,tmax/3600.], $
position=[0.15,0.30,0.85,0.50],xstyle=1,/normal

plot,hour(0:ii),alpha(0:ii),title=title,ytitle='ALPHA', $
xtitle='Time in decimal Hours',xrange=[tmin/3600.,tmax/3600.], $
position=[0.15,0.30,0.85,0.50],xstyle=1,/normal,yrange=[-1,100]

plot,hour(0:ii),psd(0:ii),title=title,ytitle='PSD at Gyrofreq', $
xtitle='Time in decimal Hours',xrange=[tmin/3600.,tmax/3600.], $
position=[0.15,0.30,0.85,0.50],xstyle=1,/normal

print, 'HOUR  ',hour(0:ii)
print, 'END of file_name loop, increment to next mass'
endfor
return
end