c	INCLUDE '../SOURCE/mod_exploit.f'

c******************************************************************

	PROGRAM des2k_diff_spl

c dessins des différences de 2 fichiers d'oscillations de type mon_modele-ad.osc
c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k

c	var: variables
c	itot: nombre de points
c	nglob: nombre de constantes global
c	nvar: nb. de variables du tableau var sans les éléments chimiques
c	nom_elem: nom des éléments du vecteur de comp. chim. généralisée
c	nchim: nombre d'éléments strictement chimiques

c	glob: variables globales du fichier mon_modele-ad.osc
c		glob(1)=mstar*msol
c		glob(2)=rtot*rsol
c		glob(3)=ltot*lsol
c		glob(4)=z0
c		glob(5)=x0
c		glob(6)=alpha
c		glob(7)=X dans ZC
c		glob(8)=Y dans ZC
c		glob(9)=d2p
c		glob(10)=d2ro
c		glob(11)=age
c		glob(12)=wrot vitesse de rotation globale
c		glob(13)=w_rot initial

c	var: variables locales utilisées

c	nvar=22 pour oscillations adiabatiques
c		var(1,i)=r*rsol
c	 	var(2,i)=log(m/mstar) -1.d38 au centre
c		var(3,i)=t
c		var(4,i)=Ptot
c		var(5,i)=ro
c		var(6,i)=gradient réel d ln T / d ln P
c		var(7,i)=l*lsol
c		var(8,i)=kap
c		var(9,i)=énergie thermo+gravifique
c		var(10,i)=grand Gamma1
c		var(11,i)=gradient adiabatique
c		var(12,i)=delta
c		var(13,i)=cp
c		var(14,i)=mu elec.
c		var(15,i)=vaissala, 0 au centre
c	 	var(16,i)=vitesse angulaire, radian/sec
c	 	var(17,i)=d ln kappa / d ln T
c	 	var(18,i)=d ln kappa / d ln ro
c	 	var(19,i)=d epsilon(nuc) / d ln T
c	 	var(20,i)=d epsilon(nuc) / d ln ro
c		var(21,i)=Ptot / Pgaz ou grad_mu
c		var(22,i)=gradient radiatif

c	  	var(22+j,i)=xchim(j)*nucleo(j), j=1,nchim

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

	USE mod_donnees, ONLY : aradia, cpturb, msol, nchim, nom_elem,
	1 nom_fich2, nom_output, nom_rot, rot_min, rsol
	USE mod_exploit, ONLY : glob, lit_nl, read_ascii, var
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, max_local
	USE mod_variables, ONLY : wrot
	
	IMPLICIT NONE

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: var1, var2
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: absc1, absc1t, absc2,
	1 absc2t, dfx1, dfx2, fx1, fx2
	REAL (kind=dp) :: pas, x

	INTEGER, PARAMETER :: ndes=1000, nf=14
	REAL (kind=sp), DIMENSION(nf,ndes) :: fdes=0.	
	REAL (kind=sp), DIMENSION(ndes) :: des=0., xdes=0.
	REAL (kind=sp), DIMENSION(nf) :: ymax=0., ymin=0.
	REAL (kind=sp) :: beta, c1, c2, xmax, xmin
	
c	INTEGER, PARAMETER :: mspl=4
	INTEGER :: i, ic, itot1, itot2, knot1, knot2, l=1, mspl, nglob, nvar,
	1 nv1, nv2
	
	LOGICAL :: en_rayon, ok
	
	CHARACTER (len=1) :: oui	
	CHARACTER (len=7) :: absc, device
	CHARACTER (len=100), DIMENSION(nf) :: titre	
	CHARACTER (len=80), DIMENSION(4) :: abid		
	CHARACTER (len=80) :: chaine, modele1, modele2
	CHARACTER (len=100) :: modelee
		
c-------------------------------------------------------------------

2000	FORMAT(8es10.3)

	PRINT*,'nom, sans extension, du premier fichier Ex: soleil'
	READ*,nom_fich2

	PRINT*,'interpolation dans le premier fichier aux points du second'
	PRINT*,'2 X (premier-second) / (premier+second)'
	
	chaine=TRIM(nom_fich2)//'.don'

	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(.NOT.ok)THEN
	 PRINT*,'Arrêt: fichier de données inconnu: ',TRIM(chaine)
	 STOP
	ENDIF
	
c lecture du fichier de données
	CALL lit_nl(wrot)
	
	PRINT*,'pour unfichier créé par:'
	PRINT*,'CESAM2k, Ex: soleil-ad.osc, entrer o'	
	PRINT*,'f037_2k, Ex: soleil_f037-ad.osc, entrer n '
	READ*,oui	
	
	IF(nom_output == 'osc_adia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-ad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-ad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_nadia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-nad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-nad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_invers')THEN	
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-inv.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-inv.osc'	
	 ENDIF
	ELSE
	 chaine=nom_output
	 WRITE(*,4)TRIM(chaine) ; STOP  
4	 FORMAT('Arrêt, pour un dessin nom_output doit être',/,
	1 'osc_adia ou osc_nadia ou osc_invers; nom_output : ',a)
	ENDIF
	
	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(ok)THEN
	 WRITE(*,5)TRIM(chaine)
5	 FORMAT('fichier ASCII utilisé: ',a)	 
	ELSE
	 WRITE(*,3)TRIM(chaine) ; STOP	
3	 FORMAT('Arrêt, fichier ASCII non trouvé: ',a)
	ENDIF
	
	CALL read_ascii(chaine,itot1,nglob,nvar,abid)
	modele1=TRIM(chaine)//' _____ ' ; modelee=TRIM(modele1)
	
c retournement du fichier var construit avec les abscisses décroissantes
	var(:,1:itot1:+1)=var(:,itot1:1:-1)
	
c transfert var --> var1
	nv1=SIZE(var,1)	
	ALLOCATE(var1(nv1,itot1)) ; var1=var	
	
c normalisation des masses et rayons
	WHERE(var1(2,:) > -1.d30)
	 var1(2,:)=EXP(var1(2,:))*glob(1)/msol
	ELSEWHERE
	 var1(2,:)=0.d0
	END WHERE
	var1(1,:)=var1(1,:)/rsol

c pour extraire des rayons, existe aussi à la fin du programme f037_2k
c	IF(.TRUE.)THEN
	IF(.FALSE.)THEN	
	 PRINT*,'extraction des rayons dans extract.R'	 	
	 OPEN(unit=35,form='formatted',status='unknown',file='extract.R')
	 WRITE(35,2002)itot1
2002	 FORMAT(i4)
	 WRITE(35,138)var1(1,:)
138	 FORMAT(5es19.12) ; CLOSE(unit=35) ; STOP
	ENDIF		

c second fichier	
	PRINT*,'nom, sans extension, du second fichier Ex: soleil'
	READ*,nom_fich2
	chaine=TRIM(nom_fich2)//'.don'
	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(.NOT.ok)THEN
	 PRINT*,'Arrêt: fichier de données inconnu: ',TRIM(chaine)
	 STOP
	ENDIF
	
c Déallocations de fichiers alloués précédemment dans lit_nl et read_ascii
	IF(ALLOCATED(nom_rot))DEALLOCATE(nom_rot,rot_min)
	DEALLOCATE(glob,nom_elem,var)
		
c lecture du fichier de données
	CALL lit_nl(wrot)
	 	 	
	PRINT*,'pour unfichier créé par:'
	PRINT*,'CESAM, Ex: soleil-ad.osc, entrer o'	
	PRINT*,'f037_2k, Ex: soleil_f037-ad.osc, entrer n '
	READ*,oui	
	
	IF(nom_output == 'osc_adia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-ad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-ad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_nadia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-nad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-nad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_invers')THEN	
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-inv.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-inv.osc'	
	 ENDIF
	ELSE
	 WRITE(*,3)TRIM(chaine) ; WRITE(*,4) ; STOP  
	ENDIF
		
	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(ok)THEN
	 WRITE(*,5)TRIM(chaine) 
	ELSE
	 WRITE(*,3)TRIM(chaine) ; STOP	
	ENDIF
	CALL read_ascii(chaine,itot2,nglob,nvar,abid)
	modele2=TRIM(chaine)//' _ _ _ '
	modelee=TRIM(modelee)//', '//TRIM(modele2)	
	
c retournement du fichier var construit avec les abscisses décroissantes
	var(:,1:itot2:+1)=var(:,itot2:1:-1)
		
c transfert var --> var2
	nv2=SIZE(var,1)	
	ALLOCATE(var2(nv2,itot2)) ; var2=var
		
c normalisation des masses et rayons
	WHERE(var2(2,:) > -1.d30)
	 var2(2,:)=EXP(var2(2,:))*glob(1)/msol
	ELSEWHERE
	 var2(2,:)=0.d0
	END WHERE
	var2(1,:)=var2(1,:)/rsol
		
c type de graphe
	PRINT*,'dessin en rayon o/n ?' ; READ(*,'(a)')oui
	en_rayon= oui == 'o'

c cas en rayon	
	IF(en_rayon)THEN
	 PRINT*,'dessin en rayon'; absc='r/Rsol'
	 ALLOCATE(absc1(itot1),absc2(itot2))
	 absc1(:)=var1(1,:) ; absc2(:)=var2(1,:)
	 PRINT*,'Zoom? o/n' ; READ(*,'(a)')oui
	 IF(oui == 'o')THEN
	  PRINT*,'entrer rmin, rmax' ; READ*,xmin,xmax
	  WRITE(*,991)xmin,xmax
991	  FORMAT('dessin entre les rayons',8es10.3)
	 ELSE
	  xmin=absc2(1) ; xmax=absc2(itot2)
	 ENDIF
	 
c cas en masse
	ELSE
	 PRINT*,'dessin en masse' ; absc='m/Msol'
	 ALLOCATE(absc1(itot1),absc2(itot2))
	 absc1(:)=var1(2,:) ; absc2(:)=var2(2,:)
	 PRINT*,'Zoom? o/n' ; READ(*,'(a)')oui
	 IF(oui == 'o')THEN
	  PRINT*,'entrer m min, m max' ; READ*,xmin,xmax
	  WRITE(*,992)xmin,xmax
992	  FORMAT('dessin entre les masses',8es10.3)
	 ELSE
	  xmin=absc2(1) ; xmax=absc2(itot2)
	 ENDIF
	ENDIF
	
c construction des B-splines
	PRINT*,"entrer l'ordre d'interpolation spline, conseillé : 4"
	READ*,mspl

	PRINT*,'interpolation B-spline ordre : ',mspl
	ALLOCATE (absc1t(itot1+mspl),absc2t(itot2+mspl),dfx1(nv1),dfx2(nv2),
	1 fx1(nv1),fx2(nv2))
	CALL bsp1dn(nv1,var1,absc1,absc1t,itot1,mspl,knot1,.FALSE.,
	1 absc1(1),l,fx1,dfx1)
	CALL bsp1dn(nv1,var2,absc2,absc2t,itot2,mspl,knot2,.FALSE.,
	1 absc2(1),l,fx2,dfx2)

c interpolation dans le fichier1 aux points du fichier2
c les abscisses des dessins
	pas=(absc2(itot2)-absc2(1))/REAL(ndes-1,dp)
	DO i=1,ndes
	 xdes(i)=absc2(1)+REAL(i-1,dp)*pas
	ENDDO	

c caractéristiques des dessins
	device='/xw'
c	WRITE(*,*)'device ? /xw, /PS, /CPS'
c	READ(5,'(a)')device
c	CALL pgbegin(0,device,2,2)
	CALL pgbegin(0,device,1,1)
	CALL pgscf(2)		!roman font
c	CALL pgsch(1.2)
c	CALL pgslw(5)

c les titres
	 titre(1)='\gD P/P'
	 titre(2)='\gD T/T'	 
	 titre(3)='\gD \gr/\gr'
	 titre(4)='\gD \gk/\gk'	 
	 titre(5)='\gD L/L'	 
	 titre(6)='\gD Vson/Vson'
	 titre(7)='\gD N2/N2'
	 titre(8)='\gD \(0583)/\(0583)'
	 titre(9)='\gD \(0583)\dad\u/\(0583)\dad'
	 titre(10)='\gD \(0583)\d\gm\u/\(0583)\d\gm'
	 titre(11)='\gD \gm\de\u/\gm\de'
	 titre(12)='\gD \gd\gm\de\u/\gd\gm\de'	 
	 titre(13)='\gD H/H'	 
	 titre(14)='\gD \gdH\d\gm\u/\gdH\d\gm\u'	 	 	 	 
	 	 	
c interpolation dans le fichier1 aux points du fichier2
	B1: DO i=1,ndes
	 x=xdes(i)	 
	 IF(x < absc1(mspl) .OR. x > absc1(itot1))CYCLE B1	 

c interpolations
	 CALL bsp1dn(nv1,var1,absc1,absc1t,itot1,mspl,knot1,.TRUE.,
	1 x,l,fx1,dfx1)
	 CALL bsp1dn(nv1,var2,absc2,absc2t,itot2,mspl,knot2,.TRUE.,
	1 x,l,fx2,dfx2)

c les écarts
	 fdes(1,i)=2.d0*(fx1(4)-fx2(4))/(fx1(4)+fx2(4))	!pression
	 fdes(2,i)=2.d0*(fx1(3)-fx2(3))/(fx1(3)+fx2(3))	!température
	 fdes(3,i)=2.d0*(fx1(5)-fx2(5))/(fx1(5)+fx2(5))	!densité
	 fdes(4,i)=2.d0*(fx1(8)-fx2(8))/(fx1(8)+fx2(8))	!opacité

c luminosité
	 IF(fx1(7)+fx2(7) /= 0.d0)THEN	 
	  fdes(5,i)=2.d0*(fx1(7)-fx2(7))/(fx1(7)+fx2(7))
	 ELSE
	  fdes(5,i)=0.d0	 
	 ENDIF
	 
c Vson
	 c1=SQRT(fx1(10)*fx1(4)/fx1(5)) ; c2=SQRT(fx2(10)*fx2(4)/fx2(5))
	 fdes(6,i)=2.d0*(c1-c2)/(c1+c2)
	 
c Vaissala
	 IF(fx1(15)+fx2(15) /= 0.d0)THEN	 
	  fdes(7,i)=2.d0*(fx1(15)-fx2(15))/(fx1(15)+fx2(15))
	 ELSE
	  fdes(7,i)=0.d0	 
	 ENDIF

	 fdes(8,i)=2.d0*(fx1(6)-fx2(6))/(fx1(6)+fx2(6))		!gradient
	 fdes(9,i)=2.d0*(fx1(11)-fx2(11))/(fx1(11)+fx2(11))	!gradient_ad
	 fdes(10,i)=2.d0*(fx1(21)-fx2(21))/(fx1(21)+fx2(21))	!gradient_µ
	 fdes(11,i)=2.d0*(fx1(14)-fx2(14))/(fx1(14)+fx2(14))	!µ_e
	 fdes(12,i)=2.d0*(dfx1(14)-dfx2(14))/(dfx1(14)+dfx2(14))!d dµ_e	 
	 fdes(13,i)=2.d0*(fx1(23)-fx2(23))/(fx1(23)+fx2(23))	!d H  
	 fdes(14,i)=2.d0*(dfx1(23)-dfx2(23))/(dfx1(23)+dfx2(23))!d dH_µ
	ENDDO B1
c	pause

c les ordonnées
	CALL max_local(fdes,nf,xdes,ndes,xmin,xmax,ymax,ymin)
	WHERE (ymax == ymin)
	 ymax=1.e-8 ; ymin=-1.e-8
	END WHERE

c les plots
	DO i=1,nf
	 CALL pgenv(xmin,xmax,ymin(i),ymax(i),0,1)
	 CALL pglabel(absc,titre(i),modelee)
	 ic=2 ; CALL pgsci(ic)
	 des(:)=fdes(i,:)
c	 CALL pgline(ndes,xdes,des)
	 CALL pgpoint(ndes,xdes,des,2)	 
	 ic=1 ; CALL pgsci(ic)	
	ENDDO

	CALL pgend

	STOP

	END PROGRAM des2k_diff_spl