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

	SUBROUTINE lum_teff_V
	
c routine PUBLIC du mod_photo

c determination de L à partir de pi, Teff, V

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

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

	USE mod_kind
	USE mod_donnees, ONLY: g, init_ctes_phys, ln10, lsol, 
	1 mbol_sun, msol, pi, sigma
	
	IMPLICIT NONE
	
	REAL (kind=dp), DIMENSION(11) :: tfe, dtfedt, dtfedfe
	REAL (kind=dp), DIMENSION(10) :: ftg, dftgdt, dftgdg
	REAL (kind=dp), SAVE ::	cte, cte1
	REAL (kind=dp) :: bc, d_bc, bmv, d_bmv,
	1 para, dpara, teff=1.d5, dteff, v, dv, l_g, dl_g,
	2 m_bol, dm_bol, l, dl, log_l, dlog_l, fe, dfe, mass, dmass,
	3 l10t, cb_p
	
	INTEGER :: tour, tourmax=10
	
	LOGICAL, SAVE :: init=.TRUE.
	LOGICAL :: lmass
	
c-------------------------------------------------------------------
	
2000	FORMAT(8es10.3)
	
	IF(init)THEN
	 init=.FALSE.
	 CALL init_ctes_phys
	 cte=5.d0/ln10
	 cte1=LOG10(4.d0*pi*sigma*g*msol/lsol)
	ENDIF

	PRINT*,' si les ecarts sont inconnus mettre 0'		
	DO WHILE(teff > 0.d0)
	 PRINT*
	 PRINT*,'entrer Teff, dTeff, sortie: 0.d0 0.d0, Ex: 5562 80'
	 READ*,teff, dteff
	 IF(teff > 0.d0)THEN
	  PRINT*,'entrer para, dpara en mas, Exemple: 82.5 0.9'
	  READ*,para,dpara
	  PRINT*,'v, dv, Exemple: 5.81 0.03'
	  READ*,v,dv
	  v=m_abs(v,para)
	  
c flower	  
	  CALL corr_bol_flower(teff,dteff,bc,d_bc,bmv,d_bmv) ; m_bol=v+bc
	  log_l=-0.4d0*(m_bol-mbol_sun) ; l=10.d0**log_l
	  PRINT*
	  PRINT*,'suivant Flower'
	  IF(dteff*dpara*dv > 0.d0)THEN
	   dv=dm_abs(dv,dpara,para); dm_bol=dv+abs(d_bc)	
	   dlog_l=-0.4d0*dm_bol ; dl=l*ln10*dlog_l	  
	   WRITE(*,"('Log L/Lsol=',es10.3,'+/-',es10.3,
	1  ', L/Lsol=',es10.3,'+/-',es10.3)")log_l,abs(dlog_l),l,abs(dl)
	   WRITE(*,"('M_bol=',es10.3,'+/-',es10.3,
	1  ', V=',es10.3,'+/-',es10.3)")m_bol,abs(dm_bol),v,abs(dv)
	  ELSE	  	  	 
	   WRITE(*,"('Log L/Lsol=',es10.3,', L/Lsol=',es10.3)")log_l,l
c	   write(*,"('M_bol=',es10.3,', V=',es10.3)")m_bol,v
	  ENDIF
	  
c lejeune	  
	  PRINT*
	  PRINT*,'Pour Lejeune, entrer [Fe/H], d[Fe/H]'
	  PRINT*,'Ex: -0.57 0.11, si inconnu, 0. 0.'
	  READ*,fe,dfe
	  CALL corr_bol_lejeune(teff,fe,tfe,dtfedt,dtfedfe)
	  mbol_sun=4.746d0	  
	  m_bol=v+tfe(10) ; log_l=-0.4d0*(m_bol-mbol_sun) ; l=10.d0**log_l
	  PRINT*,'suivant Lejeune'
	  IF(dteff*dpara*dv > 0.d0)THEN
	   dm_bol=SQRT(dv**2+(dtfedt(10)*dteff)**2+(dtfedfe(10)*dfe)**2)
	   dlog_l=-0.4d0*dm_bol ; dl=l*ln10*dlog_l	  
	   WRITE(*,"('Log L/Lsol=',es10.3,'+/-',es10.3,
	1  ', L/Lsol=',es10.3,'+/-',es10.3)")log_l,abs(dlog_l),l,abs(dl)
	   WRITE(*,"('M_bol=',es10.3,'+/-',es10.3,
	1  ', V=',es10.3,'+/-',es10.3)")m_bol,abs(dm_bol),v,abs(dv)
	  ELSE	  	  	 
	   WRITE(*,"('Log L/Lsol=',es10.3,', L/Lsol=',es10.3)")log_l,l
c	   WRITE(*,"('M_BOL=',ES10.3,', V=',ES10.3)")M_BOL,V
	  ENDIF
	  PRINT*

c bessell	  
	  PRINT*,'pour Bessell, entrer log g, dlog g ou -m/Msol et -dm/Msol'
	  PRINT*,' Ex: log g=4.44 dlog g=0.01, ou -m/Msol=-1. -dm/Msol=-0.01'
	  PRINT*,'Ex: 4.44 0.01, -1. -0.1, si inconnu entrer 100. 100.'
	  READ*,l_g,dl_g
	  IF(l_g < 0.d0)THEN
	   mass=-l_g ; dmass=-dl_g ; l_g=4.44d0 ; lmass=.true.
	  ELSEIF(l_g > 50.d0)THEN
	   l_g=4.44d0 ; dl_g=0.d0 ; lmass=.false.
	   PRINT*,'log g et masse inconnus, on prend log g=4.44 (soleil)'
	  ELSE
	   lmass=.FALSE.
	  ENDIF
	  
	  CALL corr_bol_bess(teff,l_g,ftg,dftgdt,dftgdg)
	  mbol_sun=4.74d0
	  m_bol=v+ftg(2) ; log_l=-0.4d0*(m_bol-mbol_sun) ; l=10.d0**log_l
c	  write(*,2000)m_bol
	  IF(lmass)THEN		!on tente de converger sur CbV
	   tour=0 ; l10t=cte1+4.d0*LOG10(teff)
	   B1: DO
	    tour=tour+1 ; cb_p=ftg(2) ; l_g=l10t+LOG10(mass/l)
c	    write(*,2000)cb_p,ftg(2),l_g
	    CALL corr_bol_bess(teff,l_g,ftg,dftgdt,dftgdg)
	    m_bol=v+ftg(2) ; log_l=-0.4d0*(m_bol-mbol_sun)
	    l=10.d0**log_l
c	    print*,ftg(2),cb_p,l_g
	    IF(ABS(ftg(2)-cb_p) < 0.01)EXIT B1
	    IF(tour > tourmax)THEN
	     WRITE(*,"('pas de convergence apres',i3,' iterations')")tour
	     PRINT*,'resultats incertains' ; EXIT B1
	    ENDIF
	   ENDDO B1
	   
c estimation de l'erreur sur log g
	   dl_g=SQRT((dmass/mass/ln10)**2+(0.4d0*dv)**2+
	1  (dteff*4.d0/ln10/teff)**2+(dteff*0.4*dftgdt(2))**2)
	   dl_g=dl_g/(1.d0-dftgdg(2))	   
	  ENDIF

	  PRINT*,'suivant Bessell'
	  IF(dteff*dpara*dv*dl_g > 0.d0)THEN
	   dm_bol=SQRT(dv**2+(dftgdt(2)*dteff)**2+(dftgdg(2)*dl_g)**2)	
	   dlog_l=-0.4d0*dm_bol ; dl=l*ln10*dlog_l	  
	   WRITE(*,"('Log L/Lsol=',es10.3,'+/-',es10.3,
	1  ', L/Lsol=',es10.3,'+/-',es10.3)")log_l,abs(dlog_l),l,abs(dl)
	   WRITE(*,"('M_bol=',es10.3,'+/-',es10.3,
	1  ', V=',es10.3,'+/-',es10.3)")m_bol,abs(dm_bol),v,abs(dv)
	  ELSE	  	  	 
	   WRITE(*,"('Log L/Lsol=',es10.3,', L/Lsol=',es10.3)")log_l,l
c	   WRITE(*,"('M_bol=',es10.3,', V=',es10.3)")m_bol,v
	  ENDIF
	  WRITE(*,"('Log g=',es10.3)")l_g  	  
	 ENDIF
	 PRINT*
	ENDDO
	 
	STOP
	 
	END SUBROUTINE lum_teff_V