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

	SUBROUTINE lum_teff_RIF
	
c routine PUBLIC du mod_photo
	
c determination desluminosite, temperature effective et ecart
c à partir de pi, R, I, [Fe/H]

c Auteur: P. Morel, Laboratoire Lagrange, O.C.A., CESAM2k

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

	USE mod_kind
		
	IMPLICIT NONE
		    	
	REAL (kind=dp), DIMENSION(2,4) :: jac
	REAL (kind=dp), DIMENSION(4) ::	x, dx	    	
	REAL (kind=dp), DIMENSION(2) :: f, f0, sig
	REAL (kind=dp),	PARAMETER :: h=1.d-5, unph=1.d0+h
	REAL (kind=dp) :: r, sr, i, si, para, spara, fesh, sfesh,
	1 loge, ra, ia, stor, stor0, dstor
	
	INTEGER	:: k, j
	
c---------------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

 	loge=LOG10(EXP(1.d0))
 	
	PRINT*
	PRINT*,'calcul de L à partir de para, R, I, [Fe/H]'
	PRINT*
	
	B1: DO	
	 PRINT*,'entrer pi, sigma pi en mas'
	 PRINT*,' Exemple: 752.d0 2.7d0, sortie: -1 0'
	 READ*,para,spara
	 IF(para <= 0.d0)EXIT B1 
	 PRINT*
	 WRITE(*,"('para=',es10.3,' +/-',es10.3)")para,spara

	 PRINT*
	 PRINT*,'entrer magnitude R, sigma R, Exemple: -0.01d0 0.005d0'
	 READ*,r,sr
	 PRINT*
	 WRITE(*,"('R=',es10.3,' +/-',es10.3)")r,sr

	 PRINT*
	 PRINT*,'entrer magnitude I, sigma I, Exemple: 0.69d0 0.01d0'
	 READ*,i,si 
	 PRINT*	
	 WRITE(*,"('I=',es10.3,' +/-',es10.3)")i,si
	
	 PRINT*	
	 PRINT*,'entrer [Fe/H], sigma [Fe/H], Exemple: 0.2d0 0.02d0'
	 READ*,fesh,sfesh 		
	 PRINT*	
	 WRITE(*,"('[Fe/H=]',es10.3,' +/-',es10.3)")fesh,sfesh

c variables muettes pour le calcul du jacobien
	 x(1)=para   ; x(2)=r   ; x(3)=i  ; x(4)=fesh
	 dx(1)=spara ; dx(2)=sr ;dx(3)=si ; dx(4)=sfesh

c magnitudes absolues
	 ra=m_abs(x(2),x(1)) ; ia=m_abs(x(3),x(1)) 

	 PRINT*
	 WRITE(*,"('Magnitudes absolues')")	
	 WRITE(*,"('I=',es10.3,', R=',es10.3)")ia,ra
	 PRINT*
			
	 CALL tl_rif(ra,ia,x(4),f0(1),f0(2))
					
	 PRINT*
C	 WRITE(*,2000)f0(1),f0(2)
	 WRITE(*,"('Teff=',es10.3,', L/Lsol=',es10.3)")(f0(k),k=1,2)
	 PRINT*
C	 PAUSE
		
	 DO k=1,4
	  stor0=x(k)
	  stor=stor0*unph
	  IF(stor == 0.d0)stor=h
	  dstor=stor-stor0
	  x(k)=stor	
	  ra=m_abs(x(2),x(1)) ; ia=m_abs(x(3),x(1))
	  CALL tl_rif(ra,ia,x(4),f(1),f(2))
	  x(k)=stor0
	  jac(:,k)=(f-f0)/dstor
	 ENDDO
	 
	 PRINT* ; PRINT*,'jacobien' ; WRITE(*,"(4es10.3)")(jac(:,j),j=1,4)

	 sig=matmul(abs(jac),dx)
	
c en Log
	 PRINT*
	 PRINT*,'en Log'
	 WRITE(*,1)LOG10(F0(1)),LOGE*SIG(1)/F0(1)
1	 FORMAT('Log Teff =',es10.3,' +/-',es10.3)
	 WRITE(*,3)log10(f0(2)),loge*sig(2)/f0(2)
3	 FORMAT('Log L/Lsol =',es10.3,' +/-',es10.3)

c en valeurs naturelles
	 PRINT*
	 PRINT*,'en valeurs naturelles'
	 WRITE(*,5)f0(1),sig(1)
5	 FORMAT('Teff =',es10.3,' +/-',es10.3)
	 WRITE(*,7)f0(2),sig(2)
7	 FORMAt('L/Lsol =',es10.3,' +/-',es10.3,/)
	ENDDO B1

	STOP
	
	END SUBROUTINE lum_teff_RIF