pro lit_evol,fich,nbx,age,ray,lum,teff,temp,s_c,del_s

car='rien'
unit=1
close,unit
openr,unit,fich
readf,unit,car
readf,unit,car
readf,unit,nbx
print,'dimension du tableau:',nbx
age=dblarr(nbx)
ray=dblarr(nbx)
lum=dblarr(nbx)
teff=dblarr(nbx)
temp=dblarr(nbx)
s_c=dblarr(nbx)
del_s=dblarr(nbx)
a=dblarr(7)
for i=0,nbx-1 do begin
 readf,unit,a
 age(i)=alog10(a(0)+1d-33)
 ray(i)=a(1)
 lum(i)=alog10(a(2))
 teff(i)=a(3)
 temp(i)=a(4)
 s_c(i)=a(5)/8.3e8
 del_s(i)=a(6)/8.3e8
endfor
close,unit
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro entropies

!p.thick=3

plot,[0],[0],xrange=[0,1e4],yrange=[0.7,1.2],ystyle=1,/nodata

!p.color=64
lit_evol,'results/hd209458b/1202k/alex/evol.des',nbx,age,ray,lum,teff,temp,s_c,del_s
oplot,age,s_c
oplot,age,s_c-del_s,linestyle=2

!p.color=255
lit_evol,'results/hd209458b/1623k/alex/evol.des',nbx,age,ray,lum,teff,temp,s_c,del_s
oplot,age,s_c
oplot,age,s_c-del_s,linestyle=2

!p.color=0

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro daynight

print,'daynight: find day & night temperature using flux conservation"
print,'Enter Tmean (eq. temp.)'
read,tmean
print,'Enter delta_T'
read,delt

tm4=tmean^4
told=tmean
tday=tmean+delt/2

while (abs(told-tmean) lt 0.01) do begin
	told=tday
	tday=(2*tm4-(tday-delt)^4)^0.25
endwhile

print,'Tday=',tday,'; Tnight=',tday-delt

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro evd
;-------------
print,'Calculation of evolution with day night temperature difference'
fday='results/hd209458b/1623k/alex/evol.des'
tday=1623.
fnight='results/hd209458b/1202k/alex/evol.des'
tnight=1202.
tmean=((tday^4+tnight^4)/2)^0.25
print,'Tday=',tday,'; Tnight=',tnight,'; Tmean=',tmean
fmean='results/hd209458b/1457k/alex/evol.des'


;---reads results for night & day sides & mean teq-----
lit_evol,fnight,nbx_n,age_n,ray_n,lum_n,teff_n,temp_n,s_n,ds_n
lit_evol,fday,nbx_d,age_d,ray_d,lum_d,teff_d,temp_d,s_d,ds_d
lit_evol,fmean,nbx_m,age_m,ray_m,lum_m,teff_m,temp_m,s_m,ds_m

;---computes evolution by mixing entropies at each time step dlt---
dlt=0.05		;Dlogt
tini=1		;in Log(Myr)
tmax=4	;in 	;in Log(Myr)
npts=1+(tmax-tini)/dlt
age=fltarr(npts)
r_d=fltarr(npts)
l_d=fltarr(npts)
t_d=fltarr(npts)
ss_d=fltarr(npts)
dss_d=fltarr(npts)
r_n=fltarr(npts)
l_n=fltarr(npts)
t_n=fltarr(npts)
ss_n=fltarr(npts)
dss_n=fltarr(npts)


aged=0.
agen=0.
i_d=0
i_n=0

;--All points---
for i=0,npts-1 do begin
 ;--new time step
 if (i eq 0) then begin
  age(i)=tini
  while (age_d(i_d+1) lt tini) do i_d=i_d+1
  while (age_n(i_n+1) lt tini) do i_n=i_n+1
  x_d=(tini-age_d(i_d))/(age_d(i_d+1)-age_d(i_d))
  x_n=(tini-age_n(i_n))/(age_n(i_n+1)-age_n(i_n))
 endif else begin
  age(i)=age(i-1)+dlt
  while ((age_d(i_d+1)-aged lt dlt)and(i_d lt nbx_d-2)) do i_d=i_d+1
  while ((age_n(i_n+1)-agen lt dlt)and(i_n lt nbx_n-2)) do i_n=i_n+1
  x_d=(dlt-age_d(i_d)+aged)/(age_d(i_d+1)-age_d(i_d))
  x_n=(dlt-age_n(i_n)+agen)/(age_n(i_n+1)-age_n(i_n))
 endelse

 ;--gets mean entropy
 sday=s_d(i_d)+x_d*(s_d(i_d+1)-s_d(i_d))
 snight=s_n(i_n)+x_n*(s_n(i_n+1)-s_n(i_n))
 smean=(sday+snight)/2.
; print,'snight,sday,smean',snight,sday,smean

 ;--interpolates all quantities
 i_d=max([0,i_d-2])	;just in case...
 while ((s_d(i_d) gt smean) and (i_d lt nbx_d-2)) do i_d=i_d+1
 i_d=i_d-1
 x=(smean-s_d(i_d))/(s_d(i_d+1)-s_d(i_d))
; print,x,s_d(i_d),smean,s_d(i_d+1)
 aged=age_d(i_d)+x*(age_d(i_d+1)-age_d(i_d))
 r_d(i)=ray_d(i_d)+x*(ray_d(i_d+1)-ray_d(i_d))
 l_d(i)=lum_d(i_d)+x*(lum_d(i_d+1)-lum_d(i_d))
 t_d(i)=temp_d(i_d)+x*(temp_d(i_d+1)-temp_d(i_d))
 ss_d(i)=s_d(i_d)+x*(s_d(i_d+1)-s_d(i_d))
 dss_d(i)=ds_d(i_d)+x*(ds_d(i_d+1)-ds_d(i_d))

 i_n=min([nbx_n-2,i_n+2])	;just in case...
 while ((s_n(i_n) lt smean)and(i_n gt 0)) do i_n=i_n-1
 x=(smean-s_n(i_n))/(s_n(i_n+1)-s_n(i_n))
; print,x,s_n(i_n),smean,s_n(i_n+1)
 agen=age_n(i_n)+x*(age_n(i_n+1)-age_n(i_n))
 r_n(i)=ray_n(i_n)+x*(ray_n(i_n+1)-ray_n(i_n))
 l_n(i)=lum_n(i_n)+x*(lum_n(i_n+1)-lum_n(i_n))
 t_n(i)=temp_n(i_n)+x*(temp_n(i_n+1)-temp_n(i_n))
 ss_n(i)=s_n(i_n)+x*(s_n(i_n+1)-s_n(i_n))
 dss_n(i)=ds_n(i_n)+x*(ds_n(i_n+1)-ds_n(i_n))

 print,'i,i_n,r_n,r_d,s_n,s_d,sm',$
	i,i_n,r_n(i),r_d(i),ss_n(i),ss_d(i),smean

endfor

;-------------
!p.thick=4

plot,[0],[0],xrange=[0,4],yrange=[0.,5.],ystyle=1,/nodata

oplot,age,r_n
oplot,age,r_d,linestyle=2

oplot,age_d,ray_d,linestyle=1,color=255
oplot,age_n,ray_n,linestyle=1,color=64
;oplot,age_m,ray_m,color=64

;!p.color=255
oplot,age,l_n,linestyle=2,color=64
oplot,age,l_d,linestyle=2,color=255

oplot,age_d,lum_d,linestyle=1,color=255
oplot,age_n,lum_n,linestyle=1,color=64
;oplot,age_m,lum_m,color=64

!p.color=0

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ps
set_plot,'ps'
device,/color,/bold,file='entropies.ps'
loadct,13
evd
device,/close
set_plot,'x'
loadct,0
end