PROGRAM CALIB C C This fortran code applies the calibrations proposed in the paper C "A calibration of Geneva photometry for B to G stars in terms of Teff, C logg and [M/H]" (A&AS,1996, , ). C C It allows to determine the physical parameters Teff and logg for C unreddned B to mid-G stars of the main sequence or just above it. C Metallicity is also given for cool stars ( T<7500K ), but for hotter C stars, [M/H] has to be known a priori and can take value comprised C betweem -1 and 1. C C C C We have defined three range of temperature: C C - For hot stars (Teff>11000 K), we use the reddening-free parameters X C and Y defined by Cramer and Maeder (1979) by C C X= 0.3788+1.3784U-1.2162B1-0.8498B2-0.1554V1+0.8450G C Y=-0.8288+0.3250U-2.3228B1+2.3363B2+0.7495V1-1.0865G C C They are sensitive to effective temperature and gravity respectively. C C C - For the intermediate range (8000K { Teff , log(g) } for C different metallicities ( +1 +0.5 +0.3 +0.2 +0.1 +0.0 C -0.1 -0.2 -0.3 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5 C -5.0) with in free format C - [M/H]k (B2-V1)i dj (log(g))ijk (Te)ijk C - The file bcmci contains numerical values C of the grids { (B2-V1) , m2 } --> { [M/H] , Teff } for C different gravities ( 2.5 3 3.5 4 4.5 5), with in C free format C - log(g)k (B2-V1)i m2j [M/H]ijk (Te)ijk C - The files xcycp10i, xcycp00i and xcycm10i C contain numerical values of the grids C {X,Y} --> {Teff,log(g)} for [M/H] = +1, 0 and -1 C respectively, with in free format C - Xi Yj (log(g))ij (theta)ij (Teff)ij C - The files ptpgp10i, ptpgp00i and ptpgm10i C contain numerical values of the grids C {pT,pG} --> {Teff,log(g)} for [M/H] = +1, 0 and -1 C respectively, with in free format C - pTi pGj [M/H] (log(g))ij (theta)ij (Teff)ij C C C the six colours of the Geneva photometric system REAL U,V,B1,B2,V1,G C C parameters for hot and intermediate stars REAL X,Y,pT,pG C variables for XYTL and PTGTL REAL tef3(3),log3(3),covar(2,2),tef2(3),log2(3) REAL dte2(3),dl2(3),extr2(3) REAL dte3(3),dl3(3),extr3(3) REAL MH,xMH,log C C declaration for calibration (B2-V1,d,m2) REAL m2min,m2pas C Minimal values for B2-V1,d et m2 PARAMETER (b2v1min=-0.175,dmin=+0.25,m2min=-0.66) C steps in the inverted grids PARAMETER (b2v1pas=+0.025,dpas=+0.05,m2pas=+0.02) C Maximal number of values B2-V1, d et m2. Nv > MAX(...) PARAMETER (Nb2v1=50,Nd=60,Nm2=80,Nv=80) C Number of values for log(g) and [M/H] PARAMETER (Nm=19,Ng=6) C Maximal deviation allowed during the iterative process PARAMETER (dtmax=20,dgmax=0.02,dmmax=0.02,Nitermax=6) PARAMETER (Nwarn=2) PARAMETER (NT1=Nd*Nb2v1*Nm, NT2=Nm2*Nb2v1*Ng) REAL Vb2v1(Nb2v1),Vd(Nd),Vm2(Nm2) REAL Vm(Nm),Vg(Ng) C Physical parameters and indices REAL teff,logg,metal,b2v1,d,m2 C Scatters for Te,log(g),[M/H] REAL st,sg,sm C Deviations in Te, [M/H] and log(g) during the iterative process REAL dt,dm,dg C Temporary values during the iterative process REAL teff1,teff2,logg1,logg2,metal1,metal2 REAL Mt(Nd,Nb2v1,Nm),Mg(Nd,Nb2v1,Nm),Mm(Nm2,Nb2v1,Ng) LOGICAL exist1(Nd,Nb2v1,Nm),exist2(Nm2,Nb2v1,Ng) CHARACTER cas*1,choix*1,data*20,titre*80,fileres*20 C Remarks NB and SSI CHARACTER*3 warn(Nwarn) REAL tr(Nv,7),ext,exg,exm DATA exist1/NT1*.FALSE./,exist2/NT2*.FALSE./ C Values of [M/H] for the grids (d,B2-V1) DATA Vm/-5.,-4.5,-4.,-3.5,-3.,-2.5,-2.,-1.5,-1,-0.5, &-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.5,1./ C Values of log(g) for the grids (m2,B2-V1) DATA Vg/2.5,3.0,3.5,4.0,4.5,5./ DATA p/0/,val/4/ C C OPEN(15,FILE='bcdci',STATUS='old') OPEN(16,FILE='bcmci',STATUS='old') C C Definition des vecteurs Vb2v1, Vd et Vm2 DO 1 i=1,Nb2v1 Vb2v1(i)=b2v1pas*(i-1)+b2v1min 1 CONTINUE DO 2 i=1,Nd Vd(i)=dpas*(i-1)+dmin 2 CONTINUE DO 3 i=1,Nm2 Vm2(i)=m2pas*(i-1)+m2min 3 CONTINUE C C C Reading of the file bcdci and bcmci C 4 READ(15,*,END=6)metal,b2v1,d,logg,teff i=NINT(1+(b2v1-b2v1min)/b2v1pas) j=NINT(1+(d-dmin)/dpas) DO 5 n=1,Nm IF (metal.EQ.Vm(n)) k=n 5 CONTINUE Mt(j,i,k)=teff Mg(j,i,k)=logg exist1(j,i,k)=.TRUE. GOTO 4 6 CONTINUE C 7 READ(16,*,END=9)logg,b2v1,m2,metal i=NINT(1+(b2v1-b2v1min)/b2v1pas) j=NINT(1+(m2-m2min)/m2pas) DO 8 n=1,Ng IF (logg.EQ.Vg(n)) k=n 8 CONTINUE Mm(j,i,k)=metal exist2(j,i,k)=.TRUE. GOTO 7 9 CONTINUE C C C Reading of the Geneva colours on screen for each star or in a file. C If you use a file, the first line of it has to contain a title. C C PRINT* WRITE(6,95),' Would you work star by star (y/n) ? ' 95 FORMAT(a,$) READ(*,200)cas IF (cas.NE.'y') THEN PRINT* PRINT*,'Name of the input file which contains on each line identification' PRINT*,'of the star, U V B1 B2 V1 G, E(B2-V1) and [M/H]=-1,0,1 (in free format)' READ(*,201)data PRINT*,'Name of the output file' READ(*,201) fileres OPEN(30,FILE=data,STATUS='old') OPEN(31,FILE=fileres,STATUS='unknown') READ(30,98) titre WRITE(31,98) titre WRITE(31,*) WRITE(31,*) ' (1) : calibration in grids (B2-V1,d) or (B2-V1,m2)' WRITE(31,*) ' (2) : calibration in grids (pT,pG)' WRITE(31,*) ' (3) : calibration in grids (X,Y)' WRITE(31,*) WRITE(31,*) ' SSI : Physical parameters are determined with only one' WRITE(31,*) ' iteration in the grids for cool stars' WRITE(31,*) ' NC : Physical paramaters do not converge for cool stars' WRITE(31,*) ' In this case, the process is stopped after 6 iterations' WRITE(31,*) WRITE(31,*) ' If there is no extrapolation in the grids (Te,logg), (Te,[M/H]), we display the value 0.0' WRITE(31,*) ' If the value of extrapolation is larger than 10, the value of physical parameters become absurd' WRITE(31,*) ' In this case, we change absurd values to 0.' WRITE(31,*) 98 FORMAT(A80) WRITE(31,*) WRITE(31,500) 500 FORMAT('Star ',9X,'Teff',18X,'Logg',18X,'[M/H]', & 13x,'Remarks ',3x,'Extrapolation',/,89x,'(Te,logg)', & x,'(Te,[M/H])') IF (p.EQ.0) then PRINT*,' ' PRINT*,'Which grids would you want to use ?' PRINT*,' ' PRINT*,'1) calibration in grids (b2v1,d,m2)' PRINT*,' ' PRINT*,'2) calibration in grids (pT,pG)' PRINT*,' ' PRINT*,'3) calibration in grids (X,Y)' PRINT*,' ' PRINT*,'4) calibration in grids choosen automatically by the code' PRINT*,' ' READ(*,*)val p=p+1 ENDIF ENDIF C 10 IF (cas.EQ.'y') THEN PRINT* PRINT*,'Could you give me the 6 colours of Geneva U,V,B1,B2,V1 and G (U=10 to finish)?' READ(*,*) U,V,B1,B2,V1,G IF (U.EQ.10) GOTO 90 PRINT*,'Colour excess E(B2-V1) ?' READ(*,*) Eb2v1 PRINT*,'Metallicity of the star ([M/H]=-1,0,1) ?' READ(*,*) MH PRINT*,' ' ELSE READ(30,*,END=90)ident,U,V,B1,B2,V1,G,Eb2v1,MH ENDIF C C B2-V1,d,m2,X,Y,pT,pG C b2v1 = B2-V1 d = U-B1-1.43*(B1-B2) m2 = B1-B2-0.457*(B2-V1) X = 0.3788+1.3764*U-1.2162*B1-0.8498*B2-0.1554*V1+0.8450*G Y = -0.8288+0.3235*U-2.3228*B1+2.3363*B2+0.7495*V1-1.0865*G pT = b2v1+0.1*X+0.0635*(Y-0.2*d)-0.006 pG = Y-0.2*d+1.1*(b2v1+0.1*X)+0.31 C C parameters are corrected by a possible interstellar extinction C b2v1 = b2v1 - Eb2v1 pT = pT - Eb2v1 pG = pG - 1.1*Eb2v1 C if (cas.EQ.'y') THEN write(*,27)b2v1,d,m2 27 format('B2-V1 = ',f6.3,3x,'d = ',f6.3,3x,'m2 = ',f6.3,/) write(*,28)pT,pG 28 format('pT = ',f6.3,3x,'pG = ',f6.3,/) write(*,29)X,Y 29 format('X = ',f6.3,3x,'Y = ',f6.3,/) PRINT*,' ' PRINT*,'Which grids would you want to use ?' PRINT*,' ' PRINT*,'1) calibration in grids (b2v1,d,m2)' PRINT*,' ' PRINT*,'2) calibration in grids (pT,pG)' PRINT*,' ' PRINT*,'3) calibration in grids (X,Y)' PRINT*,' ' PRINT*,'4) calibration in grids choosen automatically by the code' PRINT*,' ' READ(*,*)val endif C C C C determination of the grids C if (val.eq.4) then test1=X-1.4 test2=Y+0.04 test3=b2v1+0.029 test4=X-1.3 test5=Y+0.06 else if (val.eq.3) then test1=-1 test2=1 test4=-1 test5=1 go to 72 else if (val.eq.2) then test1=1 test2=-1 test3=-1 test4=1 test5=-1 go to 81 else if (val.eq.1) then test1=1 test2=-1 test3=1 test4=1 test5=-1 goto 82 end if C calibration (X,Y) C 72 IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) THEN if(MH.eq.0.or.MH.eq.-1.or.MH.eq.1) then CALL XYTL(X,Y,MH,thef,log,dt,dl,covar,extr) extr3(3)=extr if (extr3(3).le.10) then tef3(3)=5040./thef dte3(3)=5040.*dt/(thef*thef) log3(3)=log dl3(3)=dl rho=0 else tef3(3)=0 dte3(3)=0 log3(3)=0 dl3(3)=0 rho=0 endif else xMH=0 do 56 i=1,2 CALL XYTL(X,Y,xMH,thef,log,dt,dl,covar,extr) extr3(i)=extr if (extr3(i).le.10) then tef3(i)=5040./thef dte3(i)=5040.*dt/(thef*thef) log3(i)=log dl3(i)=dl else tef3(i)=0 dte3(i)=0 log3(i)=0 dl3(i)=0 endif if (MH.ge.0) then xMH=1 else xMH=-1 endif 56 enddo if (extr3(1).le.10.and.extr3(2).le.10) then tef3(3)= tef3(1)*(1-abs(MH))+tef3(2)*abs(MH) dte3(3)=dte3(1)*(1-abs(MH))+dte3(2)*abs(MH) log3(3)= log3(1)*(1-abs(MH))+log3(2)*abs(MH) dl3(3)=dl3(1)*(1-abs(MH))+dl3(2)*abs(MH) extr3(3)=extr3(1)*(1-abs(MH))+extr3(2)*abs(MH) else tef3(3)=0 dte3(3)=0 log3(3)=0 dl3(3)=0 endif endif goto 99 ENDIF C C calibration (pT,pG) C 81 IF (test3.le.0) THEN if(MH.eq.0.or.MH.eq.-1.or.MH.eq.1) then CALL PTGTL(pT,pG,MH,thef,log,dt,dl,covar,extr) extr2(3)=extr if (extr2(3).le.10) then tef2(3)=5040./thef dte2(3)=5040.*dt/(thef*thef) log2(3)=log dl2(3)=dl rho=0 else tef2(3)=0 dte2(3)=0 log2(3)=0 dl2(3)=0 rho=0 endif else xMH=0 do 55 i=1,2 CALL PTGTL(pT,pG,xMH,thef,log,dt,dl,covar,extr) extr2(i)=extr if (extr2(i).le.10) then tef2(i)=5040./thef dte2(i)=5040.*dt/(thef*thef) log2(i)=log dl2(i)=dl else tef2(i)=0 dte2(i)=0 log2(i)=0 dl2(i)=0 endif if (MH.ge.0) then xMH=1 else xMH=-1 endif 55 enddo if (extr2(1).le.10.and.extr2(2).le.10) then tef2(3)= tef2(1)*(1-abs(MH))+tef2(2)*abs(MH) dte2(3)=dte2(1)*(1-abs(MH))+dte2(2)*abs(MH) log2(3)= log2(1)*(1-abs(MH))+log2(2)*abs(MH) dl2(3)=dl2(1)*(1-abs(MH))+dl2(2)*abs(MH) extr2(3)=extr2(1)*(1-abs(MH))+extr2(2)*abs(MH) else tef2(3)=0 dte2(3)=0 log2(3)=0 dl2(3)=0 endif endif goto 99 ENDIF C C calibration (B2-V1,d.m2) C C First iteration : we suppose that [M/H]=0 for the determination of the C temperature effective and gravity and that log(g)=4 for the determination C of [M/H] C 82 IF (test3.gt.0) THEN C C Stars cooler than B2-V1>0.5 are not calibrate C IF (b2v1.gt.0.5) THEN IF (cas.EQ.'y') THEN PRINT*,'This star is out of the range of validity :(B2-V1) > 0.5 ' PRINT* ELSE WRITE(31,87)ident 87 FORMAT(I8,2X,'This star is out of the range of validity :(B2-V1) > 0.5') ENDIF GOTO 10 ENDIF DO 18 i=1,Nwarn warn(i)=' ' 18 CONTINUE iter=0 C Teff CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,14),exist1(1,1,14),Nd,Nb2v1, & tr,Nv,teff,dfdy,dfdx,ext) C Sigma(teff) st = SIGMA(3,dfdx,dfdy,0,0,0) C log(g) CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,14),exist1(1,1,14),Nd,Nb2v1, & tr,Nv,logg,dfdy,dfdx,exg) C Sigma(log(g)) sg = SIGMA(3,dfdx,dfdy,0,0,0) C [M/H] CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,4),exist2(1,1,4),Nm2, & Nb2v1,tr,Nv,metal,dfdy,dfdx,exm) C Sigma([M/H]) sm = SIGMA(2,dfdx,dfdy,0,0,0) C C We stop here if the uncertainty of [M/H] is too large, because it can bend C the other parameters via the iterative process C IF (b2v1.le.0.025) THEN warn(1)='SSI' GOTO 15 ENDIF C C additional iteration C DO 16 iter=1,100 dt=teff dg=logg dm=metal C C Determination of the grids with metallicity between the initial value of [M/H] C IF (metal.LT.Vm(1)) THEN imetal=1 GOTO 12 ELSE IF (metal.GT.Vm(Nm)) THEN imetal=Nm-1 GOTO 12 ENDIF DO 11 i=1,Nm-1 IF (metal.GE.Vm(i).AND.metal.LE.Vm(i+1)) THEN imetal=i GOTO 12 ENDIF 11 CONTINUE C 12 fracm=(metal-Vm(imetal))/(Vm(imetal+1)-Vm(imetal)) C CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,imetal), & exist1(1,1,imetal),Nd,Nb2v1,tr,Nv,teff1,dfdy1,dfdx1,ext) CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,imetal+1), & exist1(1,1,imetal+1),Nd,Nb2v1,tr,Nv,teff2,dfdy2,dfdx2,ext) teff = teff1 + fracm*(teff2-teff1) st = SIGMA(1,dfdx1,dfdy1,dfdx2,dfdy2,fracm) C CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,imetal), & exist1(1,1,imetal),Nd,Nb2v1,tr,Nv,logg1,dfdy1,dfdx1,exg) CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,imetal+1), & exist1(1,1,imetal+1),Nd,Nb2v1,tr,Nv,logg2,dfdy2,dfdx2,exg) logg = logg1 + fracm*(logg2-logg1) sg = SIGMA(1,dfdx1,dfdy1,dfdx2,dfdy2,fracm) C C Determination of the grids with gravity between the initial value of logg C IF (logg.LT.Vg(1)) THEN ilogg=1 GOTO 14 ELSE IF (logg.GT.Vg(Ng)) THEN ilogg=Ng-1 GOTO 14 ENDIF DO 13 i=1,Ng-1 IF (logg.GE.Vg(i).AND.logg.LE.Vg(i+1)) THEN ilogg=i GOTO 14 ENDIF 13 CONTINUE C 14 fracg=(logg-Vg(ilogg))/(Vg(ilogg+1)-Vg(ilogg)) C Determination de [M/H], ainsi que du sigma par interpolation entre C les deux grilles (m2,B2-V1) de log(g) differents C CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,ilogg), & exist2(1,1,ilogg),Nm2,Nb2v1,tr,Nv,metal1,dfdy1,dfdx1,exm) CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,ilogg+1), & exist2(1,1,ilogg+1),Nm2,Nb2v1,tr,Nv,metal2,dfdy2,dfdx2,exm) metal=metal1+fracg*(metal2-metal1) sm = SIGMA(2,dfdx1,dfdy1,dfdx2,dfdy2,fracg) C C If the differences between the last two iterations is less than 20 K for C temperature, 0.02 dex for gravity and metallicity, we stop the process here C dt=ABS(dt-teff) dg=ABS(dg-logg) dm=ABS(dm-metal) IF (dt.LE.dtmax.AND.dg.LE.dgmax.AND.dm.LE.dmmax) GOTO 17 C Do not converge after Nitermax iterations IF (iter.GE.Nitermax) THEN IF (cas.EQ.'y') THEN IF (pT-0.1.gt.0) then WRITE(*,103)dm,dt,dg 103 FORMAT(' NC dm=',F4.2,' dt=',F5.0,' dg=',F4.2) PRINT*,' On continue ? (O/N)' READ(*,200)choix ENDIF IF (choix.NE.'o'.AND.choix.NE.'O') GOTO 17 ELSE warn(2)='NC' GOTO 17 ENDIF ENDIF C 16 CONTINUE 17 CONTINUE 15 CONTINUE C C if the value of the extrapolation (ext in the grid d vs B2-V1 and exm in C the grid m2 vs B2-V1) is larger than 10, we obtain absurd values for Teff, C logg and [M/H]. For this reason, I prefer changing absurd values to 0. C IF (ext.gt.10) THEN teff=0. st=0. logg=0. sg=0. ENDIF IF (exm.gt.10) THEN metal=0. sm=0. ENDIF C C C ENDIF C 99 CONTINUE C C Results C IF (cas.EQ.'y') THEN IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) then PRINT*,'temperature and gravity are estimated' PRINT*,'by X and Y' WRITE(*,75) ifix(tef3(3)),ifix(dte3(3)),log3(3),dl3(3),extr3(3) 75 FORMAT(/' Teff = ',I6,' +-'I5,/, & ' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(X,Y)',F6.1,/) ELSE IF (test3.lt.0) THEN PRINT*,'temperature and gravity are estimated' PRINT*,'by pT and pG' WRITE(*,73)ifix(tef2(3)),ifix(dte2(3)),log2(3),dl2(3),extr2(3) 73 FORMAT(/' Teff = ',I6,' +-'I5,/, & ' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(pT,pG)',F6.1,/) ELSE IF (test3.ge.0) THEN PRINT*,'temperature, gravity and metallicity are estimated' PRINT*,'by B2-V1, d and m2' WRITE(*,71) ifix(teff),ifix(st),metal,sm,logg,sg,ext,exm,warn 71 FORMAT(/' Teff = ',I6,' +-',I5,/,' [M/H] = ',F6.2,' +-',F5.2 & ,/,' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(B2V1,d)',F7.1,/, & ' extrap(B2V1,m2)',F7.1,/,2A4,/) ENDIF ELSE IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) THEN WRITE(31,80)ident,ifix(tef3(3)),ifix(dte3(3)),log3(3),dl3(3),MH,extr3(3) 80 FORMAT(I8,5X,I6,' +-'I5,'(3)',5X, & F6.2,' +-',F5.2,'(3)',5X,F6.2,23x,F6.1) ELSE IF (test3.lt.0) THEN WRITE(31,78)ident,ifix(tef2(3)),ifix(dte2(3)),log2(3),dl2(3),MH,extr2(3) 78 FORMAT(I8,5X,I6,' +-',I5,'(2)',5X, & F6.2,' +-',F5.2,'(2)',5X,F6.2,23x,F6.1) C ELSE IF (test3.ge.0.) THEN WRITE(31,76) ident,ifix(teff),ifix(st),logg,sg, & metal,sm,warn,ext,exm 76 FORMAT(I8,5X,I6,' +-',I5,'(1)',5X,F6.2,' +-', & F5.2,'(1)',5X,F6.2,' +-',F5.2,'(1)',2x,2A4,2x,f6.1,4x,f6.1) 83 FORMAT(7f11.3) C ENDIF ENDIF C GOTO 10 C 200 FORMAT(A1) 201 FORMAT(A20) 90 STOP END C C*************************************************************************** C FUNCTION SIGMA(I,dfdx1,dfdy1,dfdx2,dfdy2,frac) C C Variance S = sigma^2 et Covariance observationnelle C sur les couleurs pour un poids P=1 (en E-6 mag^2). C (corrigee le 4.3.1993) PARAMETER (Sb2v1=40.0,Sd=186.6,Sm2=60.8) PARAMETER (Sdb2v1=18.6,Sm2b2v1=-31.3) C C Cas de Te et log(g) IF (I.EQ.1) THEN s1=SQRT(Sb2v1*dfdx1**2+Sd*dfdy1**2+2.*Sdb2v1*dfdx1*dfdy1) s2=SQRT(Sb2v1*dfdx2**2+Sd*dfdy2**2+2.*Sdb2v1*dfdx2*dfdy2) SIGMA=(s1 + frac*(s2-s1))/1000.0 C Cas de [M/H] ELSE IF (I.EQ.2) THEN s1=SQRT(Sb2v1*dfdx1**2+Sm2*dfdy1**2+2.*Sm2b2v1*dfdx1*dfdy1) s2=SQRT(Sb2v1*dfdx2**2+Sm2*dfdy2**2+2.*Sm2b2v1*dfdx2*dfdy2) SIGMA=(s1 + frac*(s2-s1))/1000.0 C Cas de la masse, ou Iter0 pour Te et logg SEULS ELSE IF (I.EQ.3) THEN s=SQRT(Sb2v1*dfdx1**2+Sd*dfdy1**2+2.*Sdb2v1*dfdx1*dfdy1) SIGMA=s/1000.0 ENDIF RETURN END C C*************************************************************************** C SUBROUTINE BISPLI(Tt,Uu, Tm,Um,Xm,Loxm,Nt,Nu, Vm,Nv & ,Xx,Dxt,Dxu, Extrap) save debut C C Routine d'interpolation pour fonction a deux variables definie C par Xm(it,iu) en certains points d'un reseau [Tm(it), Um(iu)]. C La variable LOGICAL Loxm(it,iu) dit si Xm(it,iu) est C significatif. Si on se trouve a l'interieur de la grille des points C significatifs, Extrap vaut 0. Si Extrap excede .5 ou 1. il faut C se mefier. C Xx est la valeur de la fonction en (Tt,Uu). C Dxt et Dxu sont des derivees partielles. C C On donne: C ========= C - Tt, Uu Les deux variables independantes. C - Tm(Nt), Um(Nu) Grille des variables independantes C - Xm(Nt,Nu) Valeurs (eventuellement bidon) de la C fonction en (Tm., Um.) C - Loxm(Nt,Nu)*1 Matrice logique indiquant la signification C des Xm correspondants. C - Nt, Nu Dimensions de Tm, Um, Xm et Loxm C - Vm(Nv,7) Matrice de travail C - Nv Dimension de Vm. Nv doit etre .GE.Nt et .GE.Nv C C On recoit: C ========== C - Xx Valeur de la fonction en (Tt, Uu) C - Dxt Derivee partielle p. r. a Tt en (Tt, Uu) C - Dxu Derivee partielle p. r. a Uu en (Tt, Uu) C - Extrap = 0. (Tt, Uu) se trouve a l'interieur du C reseau de (Tm., Um.) significatifs. C Sinon nombre de mailles d'extrapolation. C REAL Tm(Nt),Um(Nu),Xm(Nt,Nu),Vm(Nv,7) LOGICAL incid LOGICAL debut, Loxm(Nt,Nu) COMMON/EXTREM/Icb,Ypb,Ich,Ysh DATA debut/.TRUE./ C IF(debut)THEN debut=.FALSE. Icb=0 Ich=0 END IF mu=0 Extrap=0. ju=NRCORD(Uu,Um,Nu,incid) DO 9 iu=1,Nu C ----------- Succession d'interpolations sur Tt ----------- mt=0 DO 8 it=1,Nt IF(loxm(it,iu))THEN mt=mt+1 Vm(mt,1)=Tm(it) Vm(mt,2)=Xm(it,iu) END IF 8 CONTINUE C ------------ Pas assez de points significatifs ----------------- IF(mt.LE.1)GO TO 9 C ------------ Extrapolation trop "violente" ------------------ extrab=(Tt-Vm(1 ,1))/(Vm(1 ,1)-Vm(2 ,1)) extrah=(Tt-Vm(mt,1))/(Vm(mt,1)-Vm(mt-1,1)) extral=0. IF(iu.EQ.ju-1.OR.iu.EQ.ju+2)extral=1. IF(iu.EQ.ju .OR.iu.EQ.ju+1)extral=2. IF(extrab.GE.extral)GO TO 9 Extrap=MAX(Extrap,extrab) IF(extrah.GE.extral)GO TO 9 Extrap=MAX(Extrap,extrah) C ------------ Cas normal ----------------- mu=mu+1 Vm(mu,3)=Um(iu) CALL SPLINF(Tt,Vm(1,1),Vm(1,2),mt, Vm(1,6),Vm(1,7) & ,Vm(mu,4),Vm(mu,5),xs) 9 CONTINUE IF(mu.EQ.0)THEN C ------------ Grille trop peu fournie -------------- Extrap=1000. Xx =0. Dxt=0. Dxu=0. ELSE IF(mu.EQ.1)THEN Extrap=100. Xx =Vm(1,4) Dxt=Vm(1,5) Dxu=0. ELSE Extrap=MAX(Extrap,(Uu-Vm(1 ,3))/(Vm(1 ,3)-Vm(2 ,3))) Extrap=MAX(Extrap,(Uu-Vm(mu,3))/(Vm(mu,3)-Vm(mu-1,3))) CALL SPLINF(Uu,Vm(1,3),Vm(1,4),mu,Vm(1,6),Vm(1,7),Xx,Dxu,xs) CALL SPLINF(Uu,Vm(1,3),Vm(1,5),mu,Vm(1,6),Vm(1,7),Dxt,xp,xs) END IF END C C*************************************************************************** C FUNCTION NRCORD(Rx,Rv,Nb,Incid) C C NRCORD sert a trouver quelle composante du vecteur ordonne Rv C est egale (alors Incid est .TRUE.) ou suit immediatement Rx. C C On donne: C ========= C - Rx Nombre a situer C - Rv Vecteur ordonne C - Nb Nombre de composantes C C On recoit: C ========== C - NRCORD C - Incid Verite de Rv(NRCORD)=Rx C Declarer LOGICAL Incid et DIMENSION Rv(nn) ou nn >= Nb. C LOGICAL Incid DIMENSION Rv(Nb) C rdif=Rv(Nb)-Rv(1) rsig=SIGN(1.,rdif) NRCORD=(Nb+1)/2 IF(rdif.EQ.0)GO TO 30 C --------- Cas ou Rx sort de l'intervalle definit par Rv(1), Rv(Nb) ------- IF((Rx-Rv(Nb))*rsig.GE.0.)THEN NRCORD=Nb ELSE IF((Rx-Rv(1))*rsig.LT.0.)THEN NRCORD=0 ELSE C ------------------- Cas ou Rv a peu de composantes ----------------------- IF(Nb.LE.30)THEN DO 19 ib=1,Nb IF((Rx-Rv(ib))*rsig)11,12,19 11 NRCORD=ib-1 GO TO 30 12 NRCORD=ib GO TO 30 19 CONTINUE C --------------------------- Grand vecteur -------------------------------- ELSE nbm=Nb-1 nit=NRCORD C 20 nita=nit nit=(nit+1)/2 IF(rsig*(Rv(NRCORD)-Rx))21,30,22 21 IF(nita.LE.1)GO TO 30 NRCORD=MIN(nbm,NRCORD+nit) GO TO 20 22 NRCORD=MAX( 1,NRCORD-nit) IF(nita.GT.1)GO TO 20 END IF END IF 30 IF(NRCORD.GE.1.AND.NRCORD.LE.Nb)THEN Incid=Rx.EQ.Rv(NRCORD) ELSE Incid=.FALSE. END IF END C C*************************************************************************** C SUBROUTINE SPLINF(Xx,Xd,Yd,Nn, Ypd,Ysd, Yy,Yp,Ys) C C Procedure d'interpolation de type bicubique a travers les C points P(i) [Xd(i),Yd(i)], i=1,Nn ou les Xd(i) forment une suite C strictement monotone. C On a une fonction Yy = F(Xx). Evidemment F(Xd(i)) = Yd(i). C F(Xx) est formee d'arcs de fonctions du 3eme deg, elle est C partout continue et 2 fois derivable. La premiere et C la seconde derivee sont partout continues. La valeur de la C premiere (resp. seconde) derivee en Xx est Yp (resp. Ys). C L'algorithme de base est le meme que SPLINE et SPLINT dans C les "NUMERICAL RECIPES". C C On donne: C --------- C - Xx Valeur de la variable independante C - Xd(Nn) Abscisses des P(i) C - Yd(Nn) Ordonnees " " C - Nn Nombre de points P(i) C C On fournit: C ----------- C - Ypd(Nn) De la memoire qui contiendra les F'(Xd(i)) C - Ysd(Nn) De la memoire qui contiendra les F"(Xd(i)) C C Si on fait plusieurs interpolations avec le meme jeu de P(i) C On gagne un temps de calcul considerable en utilisant SPLINU C des le deuxieme appel. SPLINF calcule Ypd et Ysd tandis que C SPLINU utilise les Xpd et Xsd obtenus par SPLINF. C C On recoit: C ---------- C - Yy Valeur de la fonction interpolee en Xx C - Yp Valeur de la premiere derivee en Xx C - Ys Valeur de la seconde derivee en Xx. C C COMMON/EXTREM/Icb,Ypsb,Ich,Ypsh C ------------------------------- C L'algorithme laisse deux degres de liberte au probleme. C En pratique il s'agira toujours de conditions aux limites portant C soit sur la premiere, soit sur la seconde derivee. C Icb (resp. Ich) indique si on veut agir sur la premiere ou C la seconde derivee en Xx(1) (resp.(Xx(Nn)). Ypsb/h donne la valeur C que l'on veut eventuellement imposer a ladite derivee. C C - Icb/h = 0 Valeur par defaut. La seconde derivee est nulle C a l'extremite consideree. C Option recommandee si on veut extrapoler. Prudence alors! C = 1 Premiere derivee nulle. C Option recommandee si la fonction est constante C au dela ou en deca de l'extremite consideree C Si Icb/h = 0 ou 1, Ipsb/h = 0. C = 2 On impose la valeur Ipsb/h a la seconde derivee C en Xx(1/Nn). C = 3 On impose la valeur Ipsb/h a la premiere derivee C en Xx(1/Nn). C REAL Xd(Nn),Yd(Nn),Ysd(Nn),Ypd(Nn) LOGICAL incid COMMON/EXTREM/Icb,Ypsb,Ich,Ypsh C IF(Nn.LE.1)GO TO 20 dxp=Xd(2)-Xd(1) IF(dxp.EQ.0.)THEN PRINT*,'SPLINF: Xd(1)=Xd(2)=',Xd(1) RETURN END IF ypp=(Yd(2)-Yd(1))/dxp IF(Nn.EQ.2)GO TO 20 C C ---------------- Calcul des y"(Xd(i)) ------------------ C IF(Icb.LE.0.OR.Icb.GT.3)Icb=0 IF(Ich.LE.0.OR.Ich.GT.3)Ich=0 IF(Icb.LE.1)Ypsb=0. IF(Ich.LE.1)Ypsh=0. C ----------- Condition initiale en Xx(1) --------------- IF(MOD(Icb,2).EQ.0)THEN C ---------- y"(Xx(1)) est definie ----------------- Ysd(1)=Ypsb Ypd(1)=0. ELSE C ---------- y'(Xx(1)) est definie ----------------- Ysd(1)=3.*(ypp-Ypsb)/dxp Ypd(1)=-.5 END IF C ---------- Recurrences "aller" ------------------- DO 11 i=2,Nn-1 dxm=dxp dxp=Xd(i+1)-Xd(i) IF(dxp*dxm.LE.0.)THEN PRINT 101,(j,Xd(j),j=i-1,i+1) 101 FORMAT(' SPLINF',3(' Xd(',I2,') =',E12.5),' mal ordonnes') RETURN END IF dxx=dxp+dxm rap=dxm/dxx ypm=ypp ypp=(Yd(i+1)-Yd(i))/dxp ys =2.*(ypp-ypm)/dxx den=2.+rap*Ypd(i-1) Ysd(i)=(3.*ys-rap*Ysd(i-1))/den Ypd(i)=(rap-1.)/den 11 CONTINUE C ----------- Condition initiale en Xx(Nn) --------------- den=2.+Ypd(Nn-1) IF(MOD(Ich,2).EQ.0)THEN C ---------- y"(Xx(Nn)) est definie ----------------- Ysd(Nn)=Ypsh Ypd(Nn)=ypp+(den*Ypsh+Ysd(Nn-1))*dxp/6. ELSE C ---------- y'(Xx(Nn)) est definie ----------------- Ypd(Nn)=Ypsh ys =2.*(Ypsh-ypp)/dxp Ysd(Nn)=((Ypsh-ypp)*6./dxp-Ysd(Nn-1))/den END IF C -------- Calcul des y"(i) par resolution d'un systeme tridiagonal C recurrence "retour" --------------------------------- DO 12 i=Nn-1,1,-1 Ysd(i)=Ysd(i)+Ypd(i)*Ysd(i+1) Ypd(i)=Ypd(i+1)-(Xd(i+1)-Xd(i))*(Ysd(i)+Ysd(i+1))/2. 12 CONTINUE C ENTRY SPLINU(Xx,Xd,Yd,Nn,Ypd,Ysd, Yy,Yp,Ys) C C ---------------- Cas avec peu de points -------------------- 20 IF(Nn.EQ.1)THEN Ys=0. Yp=0. Yy=Yd(1) ELSE IF(Nn.EQ.2)THEN Ys=0. Yp=(Yd(2)-Yd(1))/(Xd(2)-Xd(1)) Yy=Yd(1)+(Xx-Xd(1))*Yp ELSE im=NRCORD(Xx,Xd,Nn,incid) im=MAX(1,MIN(Nn-1,im)) in=im+1 dx=Xd(in)-Xd(im) tim=(Xx-Xd(im))/dx tin=tim-1. ysm=Ysd(im) ysn=Ysd(in) Ys=ysn *tim-ysm *tin Yp=Ypd(in)*tim-Ypd(im)*tin+(ysn-ysm )*tim*tin*dx /2. Yy=Yd (in)*tim-Yd (im)*tin+(ysn+ysm+Ys)*tim*tin*dx*dx/6. END IF END C C**************************************************************************** C SUBROUTINE PTGTL(Xx,Yy,Ic, Thef,Logg, Dthef,Dlogg,Covar,Extrap) C C Calcul de Theta effective (Thef = 5040./tef1) et de Log(g) C a partir de pT et pG, parametres pour les etoiles de temperature C intermediaire.[Fe/H] est suppose connu, egal a 0 (solaire), -1 ou +1. C La routine donne une estimation des erreurs a partir de la C matrice de covariance sur les couleurs de Geneve. Ces erreurs C s'acroissent rapidement si on sort de la region ou la fonction C F: (Thef, Log(g) |----> (pT, pG) n'est plus bijective. C Dans ladite region Extrap est nulle. Ailleurs l'extrapolation C presente des risques, du moins si Extrap depasse 0.5 a 1. C C On donne : Xx, Yy Parametres de T intermediaire C Ic +1, 0, -1 [Fe/H] C ========== C On recoit: Thef, Logg 5040./Teff et log(g) C ========== Dthef, Dlogg erreurs C Covar(2,2) Matrice de covariance des erreurs C Extrap 0. si on est dans la grille de definition C nombre de mailles d'extrap. sinon. C save debut PARAMETER (Nx=60,Ny=50,Nxy=Nx*Ny,Nv=60) REAL xm(Nx),ym(Ny),Logg,them(Ny,Nx),logm(Ny,Nx),vm(Nv,7) &,Covar(2,2) LOGICAL debut,lotl(Ny,Nx) DATA debut/.TRUE./,lua/0/,lum1,lu0,lup1/10,11,12/ &,them/Nxy*0./,logm/Nxy*0./ C IF(debut)THEN debut=.FALSE. OPEN(lum1,FILE='ptpgm10i',STATUS='OLD') OPEN(lu0 ,FILE='ptpgp00i',STATUS='OLD') OPEN(lup1,FILE='ptpgp10i',STATUS='OLD') ENDIF C DO 1 ix=1,Nx DO 2 iy=1,Ny IF (Ic.EQ.0.) THEN xm(ix)=0.02*(ix-10) ym(iy)=0.02*(iy-6) ELSE IF (Ic.LT.0) THEN xm(ix)=0.02*(ix-9) ym(iy)=0.02*(iy-6) ELSE xm(ix)=0.02*(ix-10) ym(iy)=0.02*(iy-11) ENDIF 2 CONTINUE 1 CONTINUE C IF(Ic.EQ.0)THEN lul=lu0 ELSE IF(Ic.LT.0)THEN lul=lum1 ELSE lul=lup1 END IF IF(lul.NE.lua)THEN lua=lul DO 3 iy=1,Ny DO 3 ix=1,Nx 3 lotl(iy,ix)=.FALSE. 4 READ(lul,*,END=9)xd,yd,xmet,Logg,Thef,xtef IF (Ic.EQ.0.) THEN ix=NINT(50.*xd+10.) iy=NINT(50.*yd+6.) ELSE IF (Ic.LT.0.) THEN ix=NINT(50.*xd+9.) iy=NINT(50.*yd+6.) ELSE ix=NINT(50.*xd+10.) iy=NINT(50.*yd+11.) ENDIF IF(ix.LT.1.OR.ix.GT.Nx)GO TO 4 IF(iy.LT.1.OR.iy.GT.Ny)GO TO 4 lotl(iy,ix)=.TRUE. them(iy,ix)=Thef logm(iy,ix)=Logg GO TO 4 9 REWIND lul END IF C CALL BISPLI(Yy,Xx, ym,xm,them,lotl,Ny,Nx, vm,Nv & ,Thef,dtdy,dtdx,Extrap) fact=1.+2.*Extrap*Extrap Covar(1,1)=1E-4*(.4116*dtdx*dtdx+1.5064*dtdx*dtdy & +2.5896*dtdy*dtdy) Dthef=fact*SQRT(Covar(1,1)) CALL BISPLI(Yy,Xx, ym,xm,logm,lotl,Ny,Nx, vm,Nv & ,Logg,dldy,dldx,Extrap) Covar(2,2)=1E-4*(.4116*dldx*dldx+1.5064*dldx*dldy & +2.5896*dldy*dldy) Covar(1,2)=1E-4*(.4116*dtdx*dldx+.7532*dtdx*dldy+.7532*dldx*dtdy &+2.5896*dtdy*dldy) Covar(2,1)=Covar(1,2) Dlogg=fact*SQRT(Covar(2,2)) END C C******************************************************************************* C SUBROUTINE XYTL(Xx,Yy,Ic, Thef,Logg, Dthef,Dlogg,Covar,Extrap) C C Calcul de Theta effective (Thef = 5040./tef1) et de Log(g) C a partir de X et Y, parametres de Cramer & Maeder (1979). C [Fe/H] est suppose connu, egal a 0 (solaire), -1 ou +1. C La routine donne une estimation des erreurs a partir de la C matrice de covariance sur les couleurs de Geneve. Ces erreurs C s'acroissent rapidement si on sort de la region ou la fonction C F: (Thef, Log(g) |----> (X, Y) n'est plus bijective. C Dans ladite region Extrap est nulle. Ailleurs l'extrapolation C presente des risques, du moins si Extrap depasse 0.5 a 1. C C On donne : Xx, Yy Parametres de Cramer & Maeder C Ic +1, 0, -1 [Fe/H] C ========== C On recoit: Thef, Logg 5040./Teff et log(g) C ========== Dthef, Dlogg erreurs C Covar(2,2) Matrice de covariance des erreurs C Extrap 0. si on est dans la grille de definition C nombre de mailles d'extrap. sinon. C save debut PARAMETER (Nx=100,Ny=150,Nxy=Nx*Ny,Nv=150) REAL xm(Nx),ym(Ny),Logg,them(Ny,Nx),logm(Ny,Nx),vm(Nv,7) &,Covar(2,2) LOGICAL debut,lotl(Ny,Nx) DATA debut/.TRUE./,lua/0/,lum1,lu0,lup1/20,21,22/ &,them/Nxy*0./,logm/Nxy*0./ C IF(debut)THEN debut=.FALSE. OPEN(lum1,FILE='xcycm10i',STATUS='OLD') OPEN(lu0 ,FILE='xcycp00i',STATUS='OLD') OPEN(lup1,FILE='xcycp10i',STATUS='OLD') END IF C DO 1 ix=1,Nx DO 2 iy=1,Ny IF (Ic.EQ.0.) THEN xm(ix)=0.05*(ix-1) ym(iy)=0.025*(iy-11) ELSE IF (Ic.LT.0) THEN xm(ix)=0.05*ix ym(iy)=0.025*(iy-9) ELSE xm(ix)=0.05*ix ym(iy)=0.025*(iy-15) ENDIF 2 CONTINUE 1 CONTINUE C IF(Ic.EQ.0)THEN lul=lu0 ELSE IF(Ic.LT.0)THEN lul=lum1 ELSE lul=lup1 END IF IF(lul.NE.lua)THEN lua=lul DO 3 iy=1,Ny DO 3 ix=1,Nx 3 lotl(iy,ix)=.FALSE. 4 READ(lul,*,END=9)xd,yd,xmet,Logg,Thef,xtef IF (Ic.EQ.0.) THEN ix=NINT(20.*xd+1.) iy=NINT(40.*yd+11.) ELSE IF (Ic.LT.0.) THEN ix=NINT(20.*xd) iy=NINT(40.*yd+9.) ELSE ix=NINT(20.*xd) iy=NINT(40.*yd+15.) ENDIF IF(ix.LT.1.OR.ix.GT.Nx)GO TO 4 IF(iy.LT.1.OR.iy.GT.Ny)GO TO 4 lotl(iy,ix)=.TRUE. them(iy,ix)=Thef logm(iy,ix)=Logg GO TO 4 9 REWIND lul END IF C CALL BISPLI(Yy,Xx, ym,xm,them,lotl,Ny,Nx, vm,Nv & ,Thef,dtdy,dtdx,Extrap) fact=1.+2.*Extrap*Extrap Covar(1,1)=1E-4*(.8586*dtdx*dtdx+.2590*dtdx*dtdy+1.9330*dtdy*dtdy) Dthef=fact*SQRT(Covar(1,1)) CALL BISPLI(Yy,Xx, ym,xm,logm,lotl,Ny,Nx, vm,Nv & ,Logg,dldy,dldx,Extrap) Covar(2,2)=1E-4*(1.7969*dldx*dldx+.0234*dldx*dldy+2.2319*dldy* &dldy) Covar(1,2)=1E-4*(1.7969*dtdx*dldx+.0117*dtdx*dldy+.0117*dldx*dtdy &+2.2319*dtdy*dldy) Covar(2,1)=Covar(1,2) Dlogg=fact*SQRT(Covar(2,2)) END