pro readplot

print,'Reads and plots files from the chapter JUPITER''S INTERIOR'
print,'by T. Guillot, D.J. Stevenson, W.B. Hubbard and D. Saumon'
print,'in "Jupiter", eds. F. Bagenal et al. 2003'
print,''
print,'Uses IDL'
print,'Please refer to http://www.obs-nice.fr/guillot/Jupbook for instructions, or email guillot@obs-nice.fr'
print,'Thanks'
print,''

print,'The resulting figures appear as a postscript file "jupiter_figs.ps"'
print,'All quantities are in cgs units except when explicitely stated otherwise'
print,''

set_plot,'ps'
device,file='jupiter_figs.ps'

close,1
depart: 
print,'Please enter 1 character:'
print,' 0 -> to exit'
print,' 1 -> to draw plots related to Jupiter''s evolution'
print,' A to F -> to draw plots related to Jupiter''s interior (models A to F)'
car='r'
cc='r'

read,car

if (car eq '0') then goto,fin

if (car eq '1') then begin
 print,'Reads file "jupiter_evolution.dat"'
 openr,1,'jupiter_evolution.dat'
 readf,1,cc
 readf,1,cc
 readf,1,n
 a=fltarr(7,n)
 readf,1,a
 close,1

 cpar=['AGE','R/RJUP','L/LJUP','T_EFF','T_10bar','S_CENTER','DELTA_S']
 print,'Enter 1 to 6, depending on the variable to be plotted as a function of time'
 for i=1,6 do print,i,'->',cpar(i)
 read,i
 plot_oi,a(0,*),a(i,*),xtitle='Age (Myr)',ytitle=cpar(i),thick=3,charsize=1.3,xrange=[0.1,3],xstyle=1
 goto,depart
endif

fic=' '
if (car eq 'A') then fic='jupiter_A'
if (car eq 'B') then fic='jupiter_B'
if (car eq 'C') then fic='jupiter_C'
if (car eq 'D') then fic='jupiter_D'
if (car eq 'E') then fic='jupiter_E'
if (car eq 'F') then fic='jupiter_F'
if (fic eq ' ') then goto,depart

ficzel=fic+'.zel'
openr,1,ficzel
print,'Reading file',ficzel
readf,1,cc
readf,1,cc
a=fltarr(6)
b=a
c=a
d=a
readf,1,a
readf,1,cc
readf,1,b
readf,1,cc
readf,1,c
readf,1,cc
readf,1,d
print,'Z (molecular region):',b(0)
print,'Core mass (Earth masses):',a(2)
print,'Mass of heavy elements in molecular region:',d(2)
print,'Mass of heavy elements in metallic region: ',d(3)
print,'Total mass of heavy elements:              ',d(5)
close,1

ficdes=fic+'.des'
openr,1,ficdes
print,'Reading file',ficdes
for i=1,12 do readf,1,cc
readf,1,n
a=fltarr(9,n)
readf,1,a
close,1
; Test on the temperatures (not calculated in the core)
for i=0,n-1 do if (a(3,i) lt 0.) then a(3,i)=max(a(3,*))
cvar=['M/Mtot','R/Rtot','logP','logT','logRho','logKappa','logGradc','logGrada','logGradr']

print,'Choose two quantities to plot a profile'
for i=0,8 do print,i,' -> ',cvar(i)
read,i,j

plot,a(i,*),a(j,*),xtitle=cvar(i),ytitle=cvar(j),thick=3,charsize=1.3
goto,depart

fin:
device,/close

end