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

	SUBROUTINE lum_teff_VBF
c routine PUBLIC du mod_photo
c determination des luminosités, températures effectives et écarts
c a partir de pi, B, V, [Fe/H]

c Auteur: P. Morel, Département J.D. Cassini, O.C.A., Observatoire de Nice

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) :: v, sv, b, sb, para, spara, fesh, sfesh, loge,
	1 va, ba, stor, stor0, dstor

	INTEGER	:: i,j

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

2000	FORMAT(8es10.3)

 	loge=log10(exp(1.d0))

	PRINT*
	PRINT*,'calcul de L à partir de para, V, B, [Fe/H]'
	PRINT*

	b1: DO
	 PRINT*,'entrer pi, sigma pi en mas, 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 V, sigma V, Exemple: -0.01d0 0.005d0'
	 READ*,v,sv
	 PRINT*
	 WRITE(*,"('V=',es10.3,' +/-',es10.3)")v,sv

	 PRINT*
	 PRINT*,'entrer magnitude B, sigma B, Exemple: 0.69d0 0.01d0'
	 READ*,b,sb
	 PRINT*
	 WRITE(*,"('B=',es10.3,' +/-',es10.3)")b,sb

	 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)=v   ; x(3)=b  ; x(4)=fesh
	 dx(1)=spara ; dx(2)=sv ;dx(3)=sb ; dx(4)=sfesh

c magnitudes absolues
	 va=m_abs(x(2),x(1)) ; ba=m_abs(x(3),x(1))

	 PRINT*
	 WRITE(*,"('Magnitudes absolues')")
	 WRITE(*,"('B=',es10.3,', V=',es10.3)")ba,va
	 PRINT*

	 CALL tl_vbf(va,ba,x(4),f0(1),f0(2))

	 PRINT*
c	 WRITE(*,2000)f0(1),f0(2)
	 WRITE(*,"('Teff=',es10.3,', L/Lsol=',es10.3)")(f0(i),i=1,2)
	 PRINT*
c	 pause

	 DO i=1,4
	  stor0=x(i)
	  stor=stor0*unph
	  if(stor .eq. 0.)stor=h
	  dstor=stor-stor0
	  x(i)=stor
	  va=m_abs(x(2),x(1)) ; ba=m_abs(x(3),x(1))
	  call tl_vbf(va,ba,x(4),f(1),f(2))
	  x(i)=stor0
	  jac(:,i)=(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(1x,'Log Teff =',es10.3,' +/-',es10.3)
	 WRITE(*,3)log10(f0(2)),loge*sig(2)/f0(2)
3	 FORMAT(1x,'Log L/Lsol =',es10.3,' +/-',es10.3)

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

	RETURN

	END SUBROUTINE lum_teff_VBF