INCLUDE '../SOURCE/mod_nuc.f'
	
c***********************************************************************

	PROGRAM plot_reac_nuc

c liste et plot des taux interpolés des réactions thermonocléaires
c permet la comparaison avec les taux direct
c l'initialisation utilise le fichier de donnees test.don

c Auteur: P. Morel, Département J.D. Cassini, O.C.A., CESAM2k

c--------------------------------------------------------------------

	USE mod_donnees, ONLY : lit_nl, ln10, nom_fich2, w_rot
	USE mod_kind
	USE mod_nuc, ONLY : niso_tot, nreac_tot, taux_nuc
	
	IMPLICIT NONE
	
	REAL (kind=dp), DIMENSION(nreac_tot) :: rt, qt, at
	REAL (kind=dp), DIMENSION(-1:niso_tot) :: nucleot, zit	
	REAL (kind=dp) :: t
	
	INTEGER, PARAMETER :: ndes=400
	REAL (kind=sp), DIMENSION(ndes,nreac_tot) :: rdes
	REAL (kind=sp), DIMENSION(ndes) :: xdes, ydes
	REAL (kind=sp) :: pas, tdeb, tfin, xloc, xmax, xmin, ymax, ymin, yloc

	INTEGER, DIMENSION(nreac_tot,2) :: izzt
	INTEGER, DIMENSION(10):: numero=0
	INTEGER ::  i, ic, idep, ifin, j, ndess, total=2
	
	CHARACTER (len=4), DIMENSION(-1:niso_tot) :: nom_elemt
	CHARACTER (len=10) :: device
	CHARACTER (len=20), DIMENSION(nreac_tot) :: nom_react, htext

c---------------------------------------------------------------
	
2000	FORMAT(8es10.3)
2001	FORMAT(4a)
2002	FORMAT('t=',es10.3,', t9=',es10.3)
2003	FORMAT(4(es10.3,10x))

c l'nitialisation des constantes nécessite la lecture d'un fichier de donnees
	nom_fich2='test' ; CALL lit_nl(w_rot)
	PRINT*,'ATTENTION: les facteurs 1/2!, ou 1/3! sont inclus dans les taux'

c divers taux de réactions
c	IF(.TRUE.)THEN
	IF(.FALSE.)THEN
	 PRINT*
	 t=0.05d9 ; idep=22 ; ifin=24			
	 CALL taux_nuc(t,total,rt,zit,izzt,qt,nom_react,nucleot,at,nom_elemt)
	 WRITE(*,2002)t,t*1.d-9	; rt=EXP(rt)
	 WRITE(*,2001)nom_react(idep:ifin) ; WRITE(*,2003)rt(idep:ifin)

	 PRINT*
	 t=1.5d9 ; idep=62 ; ifin=63			
	 CALL taux_nuc(t,total,rt,zit,izzt,qt,nom_react,nucleot,at,nom_elemt)
	 WRITE(*,2002)t,t*1.d-9	; rt=EXP(rt)
c	 WRITE(*,2001)nom_react(idep:ifin) ; WRITE(*,2003)rt(idep:ifin)
	 WRITE(*,2001)nom_react(22),nom_react(29),nom_react(62),nom_react(63)	
	 WRITE(*,2003)rt(22),rt(29),rt(62),rt(63)
	 t=0.5d9
	 CALL taux_nuc(t,total,rt,zit,izzt,qt,nom_react,nucleot,at,nom_elemt)
	 WRITE(*,2002)t,t*1.d-9	; rt=EXP(rt)	
	 WRITE(*,2003)rt(22),rt(29),rt(62),rt(63)
	 t=0.01d9
	 CALL taux_nuc(t,total,rt,zit,izzt,qt,nom_react,nucleot,at,nom_elemt)
	 WRITE(*,2002)t,t*1.d-9	; rt=EXP(rt)	
	 WRITE(*,2003)rt(22),rt(29),rt(62),rt(63)
	ENDIF
	
c plots	
	tdeb=0.5d6 ; tfin=5.d9 ; pas=(tfin-tdeb)/(ndes-1)
	DO i=1,ndes
	 xdes(i)=tdeb+pas*(i-1)
	 t=xdes(i)
	 CALL taux_nuc(t,total,rt,zit,izzt,qt,nom_react,nucleot,at,nom_elemt)
	 rdes(i,:)=rt(:)	 
	ENDDO	
	
	xdes=LOG10(xdes*1.e-9) ; ydes(:)=rdes(:,1)
c	WRITE(*,2000)xdes ; WRITE(*,2000)ydes
	
	xmin=xdes(1) ; xmax=xdes(ndes) ; ymax=10. ; ymin=-50.
	ymax=15. ; ymin=-110.
c	xloc=-2.8 ; yloc=6.

	device='/xw'
c	WRITE(6,*)'device ? /xw, /PS, /VPS, /CPS'
c	read(5,'(a)')device
	CALL pgbegin(0,device,1,1)
	CALL pgscf(2)		!roman font
	CALL pgsch(1.2)
	CALL pgslw(2)		!epaisseur du trait
c	CALL pgscr(4,.3,.5,1.)	!amenagement du bleu

c indices (<=10)des réactions à dessiner
c	 numero(1)=64 ; numero(2)=29 ; numero(3)=17 ; numero(4)=47

	DO j=0,6
	 ic=1 ; CALL pgsci(ic) ; xloc=-2.8
	 CALL pgenv(xmin,xmax,ymin,ymax,0,0)	
	 CALL pglabel('Log\d10\u T9','Log\d10\u taux','reactions')
	 DO i=1,10
	  numero(i)=MIN(10*j+i,nreac_tot)
	 ENDDO	
	 ic=1 ; yloc=8.
	 B1: DO i=1,nreac_tot
	  IF(COUNT(i == numero) == 0)CYCLE B1
	  ydes(:)=rdes(:,i)/ln10 ; ic=ic+1 ; CALL pgsci(ic)
	  CALL pgline(ndes,xdes,ydes)
	  WRITE(*,2005)i,TRIM(nom_react(i)),qt(i)
2005	  FORMAT(i3,' : ',a,', q=',es10.3)	 
	  WRITE(htext,2004)i,TRIM(nom_react(i))
2004	  FORMAT(i3,1x,a)	 
	  yloc=yloc-5. ; CALL pgtext(xloc,yloc,htext)	 	 	 
	 ENDDO B1	
	ENDDO
		
	CALL pgend	
				
	STOP	 	 	
	 	
	END PROGRAM plot_reac_nuc