pro lit_zel,fich
common compt,ic
common param,ides,psy
common results,npp,zzx,zzy,nxp,zxx,zxy
common glob,zmol,zmet,mcore,mztot,autres
common planet,ipla
masse=[317.83,95.147]

zsol=0.0192
car='rien'
close,1
openr,1,fich+'.zel'
readf,1,car
;print,car
readf,1,car
readf,1,yz,dyz,mcore,pice,mmol,mmet
readf,1,car
readf,1,zmol,zmolmin,zmolmax,zmet,zmetmin,zmetmax
zmol=zmol/zsol
zmolmin=zmolmin/zsol
zmolmax=zmolmax/zsol
zmet=zmet/zsol
zmetmin=zmetmin/zsol
zmetmax=zmetmax/zsol
readf,1,car
readf,1,zenv,zenvmin,zenvmax,ztot,ztotmin,ztotmax
zenv=zenv/zsol
zenvmin=zenvmin/zsol
zenvmax=zenvmax/zsol
ztot=ztot/zsol
mztot=ztot*zsol*masse(ipla)
ztotmin=ztotmin/zsol
ztotmax=ztotmax/zsol
;print,mcore,zmolmin,zmolmax,zmetmin,zmetmax

if (ides eq 1) then begin
; IDES=1: On dessine (Zmol,Zmet)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 for i=1,2 do readf,1,car
 a=dblarr(10)
 zx=dblarr(27)
 zy=dblarr(27)
 ii=0
 for j=0,2 do begin
  for i=1,3 do readf,1,car
  for i=0,8 do begin
   readf,1,a
   if ((a(4) gt 0.) and (a(7) gt 0.)) then begin
    zy(ii)=a(4)/zsol
    zx(ii)=a(7)/zsol
    ii=ii+1
   endif
   if ((j eq 0) and (i eq 1)) then autres(0:1)=[a(4)/zsol,a(7)/zsol] ;Ymol-sigma
   if ((j eq 0) and (i eq 2)) then autres(2:3)=[a(4)/zsol,a(7)/zsol] ;Ymol+sigma
   if ((j eq 0) and (i eq 3)) then autres(4:5)=[a(4)/zsol,a(7)/zsol] ;YN-sigma
   if ((j eq 0) and (i eq 6)) then autres(6:7)=[a(4)/zsol,a(7)/zsol] ;YN+sigma
  endfor
  if (j eq 0) then ij=ii
 endfor
 close,1
 if (ij gt 0) then begin
  entoure,ij,zx,zy,np,xp,yp
  zxx(nxp:nxp+np-1)=xp
  zxy(nxp:nxp+np-1)=yp
  nxp=nxp+np

  entoure,ii,zx,zy,np,xp,yp
  zzx(npp:npp+np-1)=xp
  zzy(npp:npp+np-1)=yp
  npp=npp+np
 endif
; zzy(npp:npp+3)=[zmolmin,zmolmin,zmolmax,zmolmax]
; zzx(npp:npp+3)=[zmetmin,zmetmax,zmetmin,zmetmax]
; npp=npp+4
 
endif else if (ides eq 4) then begin
; IDES=4: On dessine (Mtot,Mcore)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  if (mcore le 0.1) then mcore=-0.01
; if (zmolmax ge 1.) and (zmetmax ge 0.) then begin
;  zenvmin=max([zenvmin,0.15])
;  mtotmin=zenvmin*zsol*masse(ipla)
;  mtotmax=zenvmax*zsol*masse(ipla)
;  zzy(npp:npp+3)=[mcore*(1-0.15),mcore*(1+0.20),mcore*(1+0.20),mcore*(1-0.15)]
;  zzx(npp:npp+3)=[mtotmin,mtotmin,mtotmax,mtotmax]+zzy(npp:npp+3)
;  npp=npp+4

 for i=1,2 do readf,1,car
 a=dblarr(10)
 mz=dblarr(27)
 ii=0
 for j=0,2 do begin
  for i=1,3 do readf,1,car
  for i=0,8 do begin
   readf,1,a
   if ((a(4) gt 0.) and (a(7) gt 0.)) then begin
    mz(ii)=a(4)*mmol+a(7)*mmet
    ii=ii+1
   endif
   if ((j eq 0) and (i eq 1)) then autres(0:1)=[mcore,mcore+a(4)*mmol+a(7)*mmet] ;Ymol-sigma
   if ((j eq 0) and (i eq 2)) then autres(2:3)=[mcore,mcore+a(4)*mmol+a(7)*mmet] ;Ymol+sigma
   if ((j eq 0) and (i eq 3)) then autres(4:5)=[mcore,mcore+a(4)*mmol+a(7)*mmet] ;YN-sigma
   if ((j eq 0) and (i eq 6)) then autres(6:7)=[mcore,mcore+a(4)*mmol+a(7)*mmet] ;YN+sigma
  endfor
  if (j eq 0) then ij=ii
 endfor
 close,1

 if (ij gt 0) then begin 
  mtotmin=min(mz(0:ij-1))
  mtotmax=max(mz(0:ij-1))
  zxy(nxp:nxp+3)=[mcore*(1-0.15),mcore*(1+0.20),mcore*(1+0.20),mcore*(1-0.15)]
  zxx(nxp:nxp+3)=[mtotmin,mtotmin,mtotmax,mtotmax]+zxy(nxp:nxp+3)
  nxp=nxp+4
 endif

 if (ii gt 0) then begin
  mtotmin=min(mz(0:ii-1))
  mtotmax=max(mz(0:ii-1))
  zzy(npp:npp+3)=[mcore*(1-0.15),mcore*(1+0.20),mcore*(1+0.20),mcore*(1-0.15)]
  zzx(npp:npp+3)=[mtotmin,mtotmin,mtotmax,mtotmax]+zzy(npp:npp+3)
  npp=npp+4
 endif
endif

return
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro zzel,iides,ifil
common axes,xrg,yrg
common compt,ic
common param,ides,psy
common results,npp,zzx,zzy,nxp,zxx,zxy
common glob,zmol,zmet,mcore,mztot,autres
common obs,zobs
common planet,ipla

zsol=0.0192
ides=iides
zzx=dblarr(1000)
zzy=dblarr(1000)
zxx=dblarr(1000)
zxy=dblarr(1000)
autres=fltarr(8)
if (ides eq 0) then begin
; IDES=0: On dessine ..tout..?
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ic=1
; xrg=[0,4]
; yrg=[-0.1,0.2]
    plot,[0],[0],xstyle=1,ystyle=1,xrange=xrg,yrange=yrg,/nodata
    oplot,[1,1],yrg
    oplot,[2,2],yrg
    oplot,[3,3],yrg
    oplot,xrg,[0,0]
endif else if (ides eq 1) then begin
; IDES=1: On dessine (Zmet,Zmol)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; xrg=[0,30]
; yrg=[0,12]
    plot,[0],[0],xstyle=1,ystyle=1,xrange=xrg,yrange=yrg,$
      /nodata
    oplot,[max([xrg(0),yrg(0)]),min([xrg(1),yrg(1)])],$
      [max([xrg(0),yrg(0)]),min([xrg(1),yrg(1)])],$
      linestyle=1,thick=3
endif else if (ides eq 2) then begin
; IDES=2: On dessine (Zmol,Mcore)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; xrg=[0,8]
; yrg=[0,15]
    plot,[0],[0],xstyle=1,ystyle=1,yticks=3,xrange=xrg,yrange=yrg,$
      ytitle='M!Dcore!N/M!D!20S!N!3',charsize=1.8,/nodata
endif else if (ides eq 3) then begin
; IDES=3: On dessine (Zmol,Ztot)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; xrg=[-2,8]
; yrg=[0,8]
    plot,[0],[0],xstyle=1,ystyle=1,xrange=xrg,yrange=yrg,$
      ytitle='Z!Dtot!N/Z!D!9n!N!3',/nodata
    polyfill,[xrg(0),xrg(0),zobs,zobs,xrg(0)],$
      [yrg(0),yrg(1),yrg(1),yrg(0),yrg(0)],$
      orientation=90,spacing=1.2
    polyfill,[xrg(0),xrg(0),zobs,zobs,xrg(0)],$
      [yrg(0),yrg(1),yrg(1),yrg(0),yrg(0)],$
      orientation=0,spacing=1.2
endif else if (ides eq 4) then begin
; IDES=4: On dessine (Mtot,Mcore)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; xrg=[0,50]
; yrg=[0,25]
    plot,[0],[0],xstyle=1,ystyle=1,yminor=5,xminor=5,$
      xrange=xrg,yrange=yrg,$
      /nodata 
endif

;---------------------------------------------------------------------
cppta=120
cpptb=140
cnoppta=79
cnopptb=99
;cmixa=160
;cmixb=180
;cppta=0.85*!d.n_colors
;cpptb=0.95*!d.n_colors
;cnoppta=0.85*!d.n_colors
;cnopptb=0.95*!d.n_colors

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;	TOUS LES MODELES DE JUPITER
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

vcol=[0,32,255,224,164,124,64,0,0,0,0]
ilist=0
car='rien'
close,2
openr,2,'liste_jup.txt'
readf,2,car
while (strlen(car) gt 1) do begin
    print,'Nom de la liste: ',car,'; couleur:',vcol(ilist)
    npp=0
    nxp=0
    fic='rien'
    readf,2,fic
    while (strlen(fic) gt 1) do begin
        print,'Fichier a lire: "',fic,'"',strlen(fic)
        lit_zel,fic
        readf,2,fic
    endwhile
    entoure,npp,zzx,zzy,np,xp,yp
    polyfill,xp,yp,$
      orientation=45,thick=2,spacing=0.5,$
      clip=[xrg(0),yrg(0),xrg(1),yrg(1)],noclip=0,color=vcol(ilist)
    oplot,xp,yp,thick=3,color=vcol(ilist)
    readf,2,car
    ilist=ilist+1
endwhile
close,2

!p.thick=1

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro pttlabel,car
; Ecrit un petit label correspondant au modele
common param,ides,psy
common glob,zmol,zmet,mcore,mztot,autres

if (ides eq 1) then begin
    xyouts,[zmet],[zmol],car,charsize=1.0,align=0.5,charthick=5
endif
if (ides eq 4) then begin
    xyouts,[mztot],[mcore],car,charsize=1.0,align=0.5,charthick=5
endif
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro essai

n=7
x=[2,0,5,3,1,3,2]
y=[1,3,2,3,1,4,2]
plot,x,y,psym=4,thick=2,xstyle=2,ystyle=2

entoure,n,x,y,np,xp,yp
oplot,xp,yp

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro entoure,n,x,y,np,xp,yp

if (n le 0) then return
; We introduce a small random scatter
for i=0,n-1 do begin
 x(i)=x(i)*(1+(randomu(t)-0.5)*1e-5)
 y(i)=y(i)*(1+(randomu(t)-0.5)*1e-5)
endfor

i=0
alpha=dblarr(n+1)
betta=dblarr(n+1)
ind=intarr(n+1)
while (x(i) gt min(x(0:n-1))) do i=i+1
imin=i
c=0
alpha(c)=!pi/2
betta(c)=!pi/2
ind(c)=-1

;print,c,imin,x(imin),y(imin)
;if (n gt 30) then oplot,[x(imin)],[y(imin)],psym=4
while (ind(c) ne imin) do begin
 for j=0,n-1 do begin
  if (j ne i) then begin
   bet=atan((y(j)-y(i))/(x(j)-x(i)+1d-30))
   alph=bet-betta(c)
   if (alph lt 0) then alph=alph+!pi
   if (alph gt alpha(c+1)) then begin
    ind(c+1)=j
    betta(c+1)=bet
    alpha(c+1)=alph
   endif
  endif
 endfor
 c=c+1
 i=ind(c)
; print,c,ind(c),x(ind(c)),y(ind(c))
;if (n gt 30) then oplot,[x(ind(c))],[y(ind(c))],psym=2
endwhile
ind(0)=imin

np=c+1
xp=dblarr(np)
yp=dblarr(np)
for i=0,np-1 do begin
 xp(i)=x(ind(i))
 yp(i)=y(ind(i))
endfor

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro abundances,yhe,sol,ztot,z

c=['H','He','C','N','O','Ne','Na','Mg','Al','Si','P','S','Ar','Ca',$
	'Cr','Fe','Ni']
m=[1.,4.,12.,14.,16.,20.,23.,24.,27.,28.,31.,32.,40.,40.,52.,56.,59.]
x=[2.79e10,2.72e9,1.01e7,3.13e6,2.38e7,3.44e6,5.74e4,1.074e6,8.49e4,1e6,$
	1.04e4,5.15e5,1.01e5,6.11e4,1.35e4,9e5,4.93e4]

if (yhe ne 0) then xhe=m(0)/m(1)/(1-yhe)*yhe else xhe=x(1)/x(0)

somximi=1+m(1)*xhe
for i=2,16 do somximi=somximi+sol(i)*m(i)*x(i)/x(0)

ztot=(somximi-1-m(1)*xhe)/somximi
for i=0,16 do z(i)=sol(i)*x(i)/x(0)*m(i)/somximi

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro plo_abun,xrg,yrg,zobs
zsol=0.0192
sol=fltarr(17)
z=fltarr(17)

;xrg=[-2,8]
;yrg=[0,4]
;plot,[0],[0],xstyle=1,ystyle=1,xrange=xrg,yrange=yrg,yticks=1,yminor=1,$
;	ytickname=[' ',' '],xtitle='Z!Dmol!N/Z!D!9n!N!3',xticklen=0.16,$
;	charsize=1.8,/nodata
; everything is 3 times solar --------------------
xx=-[0,0,1,1,0]*0.5+xrg(1)
sol(*)=3.
abundances,0.,sol,ztot,z
ztot=ztot/zsol
z=z/zsol
print,ztot,z
polyfill,xx,[0,ztot,ztot,0,0],col=223
oplot,xx,[0,ztot,ztot,0,0]
xyouts,7.35,(ztot+z(2)+z(3)+z(4))/2-.2,'others',charsize=1.4,charthick=3,$
	align=1.0,col=223
polyfill,xx,[0,z(2)+z(3)+z(4),z(2)+z(3)+z(4),0,0],col=0.5*!d.n_colors
oplot,xx,[0,z(2)+z(3)+z(4),z(2)+z(3)+z(4),0,0]
xyouts,7.35,z(2)+z(3)+z(4)/2-.2,'O',charsize=1.5,charthick=3,$
	align=1.0,col=0.5*!d.n_colors
polyfill,xx,[0,z(2)+z(3),z(2)+z(3),0,0]
oplot,xx,[0,z(2)+z(3),z(2)+z(3),0,0]
xyouts,7.00,z(2)+z(3)/2-.2,'N',charsize=1.5,charthick=3,align=1.0
polyfill,xx,[0,z(2),z(2),0,0],col=111
oplot,xx,[0,z(2),z(2),0,0]
xyouts,7.35,z(2)/2-.2,'C',charsize=1.5,charthick=3,align=1.0,col=111

zobs=ztot
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro arrow,xrg,yrg,ra,xx,yy,rl

print,xx(0),yy(0),xx(1),yy(1),format='("Fleche: (",f8.5,",",f8.5,") -> (",f8.5,",",f8.5,")")'
oplot,xx,yy
xstar=ra*(xx-xrg(0))/(xrg(1)-xrg(0))
ystar=(yy-yrg(0))/(yrg(1)-yrg(0))
a=atan((ystar(1)-ystar(0)+1e-33)/(xstar(1)-xstar(0)+1e-33))
if (xx(1) lt xx(0)) then a=a+!pi
xa=fltarr(3)
ya=fltarr(3)
xa(0)=xstar(1)-rl*sin(!pi/4+a)
ya(0)=ystar(1)+rl*cos(!pi/4+a)
xa(1)=xstar(1)
ya(1)=ystar(1)
xa(2)=xstar(1)-rl*sin(!pi/4-a)
ya(2)=ystar(1)-rl*cos(!pi/4-a)
xa=xrg(0)+xa*(xrg(1)-xrg(0))/ra
ya=yrg(0)+ya*(yrg(1)-yrg(0))
oplot,xa,ya,linestyle=0
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro pszel
fic='zeljupsat.ps'
set_plot,'ps'
device,ysize=23,yoffset=3,xsize=17,xoffset=2.0,scale=1,/color,file=fic
loadct,13
zel
device,/close
set_plot,'x'
print,'J''ai cree le fichier "',fic,'"'
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro zel
; JUPITER ONLY ;
common axes,xrg,yrg
common obs,zobs
common planet,ipla
loadct,13

ipla=0
!p.charthick=2
erase
xrg=[0,7]
yrg=[0,7]
!x.ticklen=0.04
!x.title='Z!Dmet!N/Z!D!9n!N!3'
!y.title='Z!Dmol!N/Z!D!9n!N!3'
!x.minor=2
!y.minor=2
!p.noerase=1
!p.charsize=1.8
!p.position=[.20,.15,.90,.49]
zzel,1,1
zzel,1,-1

xrg=[0,50]
yrg=[0,15]
!x.ticklen=0.04
!p.charsize=1.8
!x.title='M!DZtot!N/M!D!20S!N!3'
!y.title='M!Dcore!N/M!D!20S!N!3'
!x.minor=0
!y.minor=0
!p.noerase=1
!p.position=[.20,.60,.90,.95]
zzel,4,1
zzel,4,0
xyouts,.87,.90,'Jupiter',charthick=3,align=1.0,/normal

end