SUBROUTINE BVRI(MODE,FE,GV,TEFF,BMV,VMR,VMI,BCV) C ---------------------------------------------------------------------- C *** THIS READS, AND INTERPOLATES IN, bvrilo.data AND bvrihi.data (SEE C VANDENBERG & CLEM 2003, AJ, 126, IN PRESS). C *** SET MODE=0 TO READ THE TABLES AND TO INTERPOLATE THEM TO THE C DESIRED VALUE OF [Fe/H] (= FE IN THE ARGUMENT LIST) C *** SET MODE=1 TO INTERPOLATE FOR THE COLORS AT THE DESIRED VALUES OF C [Fe/H], log g, and log Teff (GV AND TEFF REPRESENT THE C LAST TWO OF THESE QUANTITIES IN THE ARGUMENT LIST). C (NOTE THAT THE [Fe/H] INTERPOLATION IS CARRIED OUT ONLY C WHEN THE SUBROUTINE IS CALLED WITH MODE=0 OR MODE=-1.) C *** SET MODE=-1 IF THE INPUT TABLES, WHICH HAVE ALREADY BEEN READ C USING MODE=0, ARE TO BE RE-INTERPOLATED TO BE C CONSISTENT WITH A NEW VALUE OF [Fe/H]. C *** THE OUTPUT CONSISTS OF THE B-V, V-R, AND V-I COLORS, ALONG WITH C THE BOLOMETRIC CORRECTION TO V (ON THE SCALE WHERE THE SUN HAS C M_bol = 4.75 and M_V = 4.84). THESE QUANTITIES ARE REPRESENTED, C IN TURN, BY BMV, VMR, VMI, AND BCV IN THE ARGUMENT LIST. (B AND V C ARE ON THE JOHNSON SYSTEM, R AND I ON THE COUSINS SYSTEM.) C *** MODE MUST BE SET TO ZERO THE FIRST TIME THAT THIS SUBROUTINE IS C CALLED, AS IT IS ONLY WHEN MODE=0 THAT BVRILO.DATA AND BVRIHI.DATA C ARE READ. ONCE THE TABLES HAVE BEEN INTERPOLATED TO A PARTICULAR C VALUE OF [Fe/H] (USING MODE=0 OR MODE=-1), INTERPOLATIONS FOR THE C COLORS APPROPRIATE TO INPUT VALUES OF [Fe/H], log g, and log Teff C MAY BE CARRIED OUT *ANY NUMBER OF TIMES* USING MODE=1 (THE ONLY C VALUE OF MODE FOR WHICH COLOR INFORMATION IS OBTAINED). C *** NOTE THAT THE BVRILO.DATA AND BVRIHI.DATA FILES ARE "ATTACHED" TO C UNITS 10 AND 11, RESPECTIVELY. WHENEVER AN [Fe/H] INTERPOLATION C IS CARRIED OUT, A MESSAGE IS SENT TO UNIT 6 ("ATTACHED" TO THE C TERMINAL) TO INDICATE THE VALUE OF [FE/H] THAT APPLIES TO C SUBSEQUENT INTERPOLATIONS FOR log g and log Teff. C *** CHKBVRI.FOR (ALSO ON DISK) MAY BE USED TO TEST THE INTERPOLATION C SUBROUTINE C *** NOTE, IN THE OPEN STATEMENTS BELOW, THAT THE NAMES OF THE INPUT C *** DATA FILES ARE WRITTEN IN LOWER-CASE TEXT C ---------------------------------------------------------------------- *23456789012345678901234567890123456789012345678901234567890123456789012 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(34,7,6,4),B(34,7,4),BV(11,12,6),VR(11,12,6), 1 VI(11,12,6),BC(11,12,6),C(11,12),D(11,12),E(11,12),F(11,12), 2 H(4,4),TB(34),TT(34),GG(7),TA(11),T(11),G(12),FEH(6),P(4),Q(4), 3 R(4),S(4),AG(4),AT(4),X(11),Y(11) CHARACTER*19 NAMCOL CHARACTER*10 NAMFE SAVE 1000 FORMAT(I3,13X,I3,13X,I3,13X,I3) 1001 FORMAT(13F6.0) 1002 FORMAT(20F4.1) 1003 FORMAT(13F6.3) 1004 FORMAT(' Color grid interpolated to [FE/H] =',F6.2) 1005 FORMAT(' log g =',F6.3,' log Teff =',F7.4,' ARE OUTSIDE THE', 1 1X,'RANGE OF THE COLOR TABLES') 1006 FORMAT(A10,F5.2,A19) 1007 FORMAT(1X,' **** INPUT DATA FILE DOES NOT EXIST **** ') IF(MODE) 11,3,16 C *** WHEN MODE=0 THE COLOR TRANSFORMATION TABLES ARE READ 3 TBND=LOG10(5000.) GBND=2. OPEN(UNIT=10,ERR=49,FILE='bvrilo.data',STATUS='OLD') OPEN(UNIT=11,ERR=49,FILE='bvrihi.data',STATUS='OLD') READ(10,1000) NT,NG,NFE,NDX READ(10,1001) (TA(I),I=1,NT) READ(10,1002) (G(I),I=1,NG) DO 4 I=1,NT 4 T(I)=LOG10(TA(I)) DO 5 K=1,NFE READ(10,1006) NAMFE,FEH(K),NAMCOL DO 5 J=1,NG READ(10,1003) (BV(I,J,K),I=1,NT) 5 CONTINUE DO 6 K=1,NFE READ(10,1006) NAMFE,FEH(K),NAMCOL DO 6 J=1,NG READ(10,1003) (VR(I,J,K),I=1,NT) 6 CONTINUE DO 7 K=1,NFE READ(10,1006) NAMFE,FEH(K),NAMCOL DO 7 J=1,NG READ(10,1003) (VI(I,J,K),I=1,NT) 7 CONTINUE DO 8 K=1,NFE READ(10,1006) NAMFE,FEH(K),NAMCOL DO 8 J=1,NG READ(10,1003) (BC(I,J,K),I=1,NT) 8 CONTINUE READ(11,1000) NTT,NGG,NFE,NDXX READ(11,1001) (TB(I),I=1,NTT) READ(11,1002) (GG(I),I=1,NGG) DO 9 I=1,NTT 9 TT(I)=LOG10(TB(I)) DO 10 L=1,NDXX DO 10 K=1,NFE READ(11,1006) NAMFE,FEH(K),NAMCOL DO 10 J=1,NGG READ(11,1003) (A(I,J,K,L),I=1,NTT) 10 CONTINUE CLOSE(UNIT=10,STATUS='KEEP') CLOSE(UNIT=11,STATUS='KEEP') C *** WHEN MODE=0 OR MODE=-1, COLOR TRANSFORMATION TABLES ARE C *** CREATED FOR THE INPUT [Fe/H] VALUE USING LINEAR INTERPOLATION 11 DO 12 M=2,NFE K=M-1 IF(FE.LE.FEH(M)) GO TO 13 12 CONTINUE 13 M=K+1 SLOPE=(FE-FEH(K))/(FEH(M)-FEH(K)) DO 14 J=1,NG DO 14 I=1,NT C(I,J)=BV(I,J,K)+SLOPE*(BV(I,J,M)-BV(I,J,K)) D(I,J)=VR(I,J,K)+SLOPE*(VR(I,J,M)-VR(I,J,K)) E(I,J)=VI(I,J,K)+SLOPE*(VI(I,J,M)-VI(I,J,K)) F(I,J)=BC(I,J,K)+SLOPE*(BC(I,J,M)-BC(I,J,K)) 14 CONTINUE DO 15 L=1,NDXX DO 15 J=1,NGG DO 15 I=1,NTT B(I,J,L)=A(I,J,K,L)+SLOPE*(A(I,J,M,L)-A(I,J,K,L)) 15 CONTINUE WRITE(6,1004) FE C *** WHEN MODE=0 OR MODE=-1, CONTROL RETURNS TO THE CALLING C *** PROGRAM ONCE THE [Fe/H] INTERPOLATION IS CARRIED OUT (I.E., NO C *** INTERPOLATIONS ARE PERFORMED FOR INPUT VALUES OF GV AND TEFF) GO TO 50 C *** WHEN MODE=1, INTERPOLATIONS ARE CARRIED OUT FOR THE INPUT VALUES C *** OF [Fe/H], log g, AND log Teff (= FE, GV, AND TEFF IN THE ARGUMENT C *** LIST) 16 IF(TEFF.LE.TBND) GO TO 18 IF(GV.GE.GBND) GO TO 40 IF(TEFF.LE.T(NT)) GO TO 18 WRITE(8,1005) GV,TEFF C *** EXECUTION HALTS WITH A "STOP30" CODE IF THE INPUT TEMPERATURE IS C *** OUTSIDE THE RANGE OF THE TABLES ON THE HIGH SIDE STOP30 C *** THE NEXT SECTION ASSUMES THAT THE LOW-TEMPERATURE TABLES ARE THE C *** RELEVANT ONES TO USE 18 NGM=NG-1 DO 19 I=3,NGM MG=I-2 IF(GV.LE.G(I)) GO TO 20 19 CONTINUE IF(GV.LE.G(NG)) GO TO 20 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log g CONSIDERED IN THE C *** TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT TWO C *** LINES SHOULD BE ACTIVATED C WRITE(8,1005) GV,TEFF C STOP31 20 R(1)=G(MG) R(2)=G(MG+1) R(3)=G(MG+2) R(4)=G(MG+3) CALL LGRAN4(R,AG,GV) NTM=NT-1 DO 25 I=3,NTM MT=I-2 IF(TEFF.LE.T(I)) GO TO 30 25 CONTINUE IF(TEFF.LE.T(NT)) GO TO 30 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log Teff CONSIDERED IN C *** THE TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT C *** TWO LINES SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP32 30 R(1)=T(MT) R(2)=T(MT+1) R(3)=T(MT+2) R(4)=T(MT+3) CALL LGRAN4(R,AT,TEFF) L=MG-1 DO 35 I=1,4 L=L+1 P(I)=AT(1)*C(MT,L)+AT(2)*C(MT+1,L)+AT(3)*C(MT+2,L)+AT(4)*C(MT+3,L) Q(I)=AT(1)*D(MT,L)+AT(2)*D(MT+1,L)+AT(3)*D(MT+2,L)+AT(4)*D(MT+3,L) R(I)=AT(1)*E(MT,L)+AT(2)*E(MT+1,L)+AT(3)*E(MT+2,L)+AT(4)*E(MT+3,L) 35 S(I)=AT(1)*F(MT,L)+AT(2)*F(MT+1,L)+AT(3)*F(MT+2,L)+AT(4)*F(MT+3,L) C *** 4-POINT LAGRANGIAN INTERPOLATION IS USED TO DERIVE THE VALUES OF C *** B-V, V-R, V-I, AND BC_V CORRESPONDING TO THE INPUT VALUES OF C *** [Fe/H], log g, and log Teff. BMV=AG(1)*P(1)+AG(2)*P(2)+AG(3)*P(3)+AG(4)*P(4) VMR=AG(1)*Q(1)+AG(2)*Q(2)+AG(3)*Q(3)+AG(4)*Q(4) VMI=AG(1)*R(1)+AG(2)*R(2)+AG(3)*R(3)+AG(4)*R(4) BCV=AG(1)*S(1)+AG(2)*S(2)+AG(3)*S(3)+AG(4)*S(4) C *** SPLINE INTERPOLATION IS USED TO FIND THE VALUE OF BC_V AT LOW C *** TEMPERATURES. DO 38 I=1,NT X(I)=T(I) 38 Y(I)=AG(1)*F(I,MG)+AG(2)*F(I,MG+1)+AG(3)*F(I,MG+2)+AG(4)*F(I,MG+3) CALL INTEP(TEFF,BCV,X,Y,NT,IER) GO TO 50 C *** THE NEXT SECTION ASSUMES THAT THE HIGH-TEMPERATURE TABLES ARE THE C *** RELEVANT ONES TO USE. 40 IF(TEFF.LE.TT(NTT)) GO TO 42 WRITE(8,1005) GV,TEFF STOP33 42 NGM=NGG-1 DO 43 I=3,NGM MG=I-2 IF(GV.LE.GG(I)) GO TO 44 43 CONTINUE IF(GV.LE.GG(NGG)) GO TO 44 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log g CONSIDERED IN THE C *** TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT TWO C *** LINES SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP34 44 R(1)=GG(MG) R(2)=GG(MG+1) R(3)=GG(MG+2) R(4)=GG(MG+3) CALL LGRAN4(R,AG,GV) NTM=NTT-1 DO 45 I=3,NTM MT=I-2 IF(TEFF.LE.TT(I)) GO TO 46 45 CONTINUE IF(TEFF.LE.TT(NTT)) GO TO 46 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log Teff CONSIDERED IN C *** THE TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT C *** TWO LINES SHOULD SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP35 46 R(1)=TT(MT) R(2)=TT(MT+1) R(3)=TT(MT+2) R(4)=TT(MT+3) CALL LGRAN4(R,AT,TEFF) L=MG-1 DO 48 I=1,4 L=L+1 DO 48 K=1,NDXX H(I,K)=AT(1)*B(MT,L,K)+AT(2)*B(MT+1,L,K)+AT(3)*B(MT+2,L,K)+ 1 AT(4)*B(MT+3,L,K) 48 CONTINUE C *** 4-POINT LAGRANGIAN INTERPOLATION IS USED TO DERIVE THE VALUES OF C *** B-V, V-R, V-I, AND BC_V CORRESPONDING TO THE INPUT VALUES OF C *** [Fe/H], log g, and log Teff. BMV=AG(1)*H(1,1)+AG(2)*H(2,1)+AG(3)*H(3,1)+AG(4)*H(4,1) VMR=AG(1)*H(1,2)+AG(2)*H(2,2)+AG(3)*H(3,2)+AG(4)*H(4,2) VMI=AG(1)*H(1,3)+AG(2)*H(2,3)+AG(3)*H(3,3)+AG(4)*H(4,3) BCV=AG(1)*H(1,4)+AG(2)*H(2,4)+AG(3)*H(3,4)+AG(4)*H(4,4) GO TO 50 49 WRITE(6,1007) STOP 50 RETURN END SUBROUTINE UVBY(MODE,FE,GV,TEFF,SBY,SM1,SC1,BCV) C ---------------------------------------------------------------------- C *** THIS READS, AND INTERPOLATES IN, UVBYLO.DATA AND UVBYHI.DATA (SEE C CLEM, VANDENBERG, GRUNDAHL, & BELL, AJ, IN PRESS). C *** SET MODE=0 TO READ THE TABLES AND TO INTERPOLATE THEM TO THE C DESIRED VALUE OF [Fe/H] (= FE IN THE ARGUMENT LIST) C *** SET MODE=1 TO INTERPOLATE FOR THE COLORS AT THE DESIRED VALUES OF C [Fe/H], log g, and log Teff (GV AND TEFF REPRESENT THE C LAST TWO OF THESE QUANTITIES IN THE ARGUMENT LIST). C (NOTE THAT THE [Fe/H] INTERPOLATION IS CARRIED OUT ONLY C WHEN THE SUBROUTINE IS CALLED WITH MODE=0 OR MODE=-1.) C *** SET MODE=-1 IF THE INPUT TABLES, WHICH HAVE ALREADY BEEN READ C USING MODE=0, ARE TO BE RE-INTERPOLATED TO BE C CONSISTENT WITH A NEW VALUE OF [Fe/H]. C *** THE OUTPUT CONSISTS OF THE b-y, m1, AND c1 INDICES, ALONG C WITH THE BOLOMETRIC CORRECTION TO V (ON THE SCALE WHERE THE SUN C HAS M_bol = 4.75 and M_V = 4.84). THESE QUANTITIES ARE C REPRESENTED, IN TURN, BY SBY, SM1, SC1, AND BCV IN THE C ARGUMENT LIST. NOTE THAT IT IS ASSUMED THE BOLOMETRIC C CORRECTIONS TO STROMGREN y ARE IDENTICAL TO THOSE IN JOHNSON V. C *** MODE MUST BE SET TO ZERO THE FIRST TIME THAT THIS SUBROUTINE IS C CALLED, AS IT IS ONLY WHEN MODE=0 THAT UVBYLO.DATA AND UVBYHI.DATA C ARE READ. ONCE THE TABLES HAVE BEEN INTERPOLATED TO A PARTICULAR C VALUE OF [Fe/H] (USING MODE=0 OR MODE=-1), INTERPOLATIONS FOR THE C COLORS APPROPRIATE TO INPUT VALUES OF [Fe/H], log g, and log Teff C MAY BE CARRIED OUT *ANY NUMBER OF TIMES* USING MODE=1 (THE ONLY C VALUE OF MODE FOR WHICH COLOR INFORMATION IS OBTAINED). C *** NOTE THAT THE UVBYLO.DATA AND UVBYHI.DATA FILES ARE "ATTACHED" TO C UNITS 12 AND 13, RESPECTIVELY. WHENEVER AN [Fe/H] INTERPOLATION C IS CARRIED OUT, A MESSAGE IS SENT TO UNIT 6 ("ATTACHED" TO THE C TERMINAL) TO INDICATE THE VALUE OF [FE/H] THAT APPLIES TO C SUBSEQUENT INTERPOLATIONS FOR log g and log Teff. C *** CHKUVBY.FOR (ALSO ON DISK) MAY BE USED TO TEST THE INTERPOLATION C SUBROUTINE C ---------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 M1(13,12,8) DIMENSION A(51,7,8,4),B(51,7,4),BY(13,12,8),C1(13,12,8), 1 BC(13,12,8),C(13,12),D(13,12),E(13,12),F(13,12),H(4,4),TB(51), 2 TT(51),GG(7),TA(13),T(13),G(12),FEH(8),O(4),P(4),Q(4),R(4), 3 AG(4),AT(4),X(13),Y(13) CHARACTER*19 NAMCOL CHARACTER*10 NAMFE SAVE 1000 FORMAT(I3,13X,I3,13X,I3,13X,I3) 1001 FORMAT(13F6.0) 1002 FORMAT(20F4.1) 1003 FORMAT(13F6.3) 1004 FORMAT(' Color grid interpolated to [FE/H] =',F6.2) 1005 FORMAT(1X,7HLOG G =,F6.3,11H LOG TEFF =,F7.4,8H OUTSIDE, 1 1X,11HCOLOR TABLE) 1006 FORMAT(A10,F5.2,A19) 1007 FORMAT(1X,' **** INPUT DATA FILE DOES NOT EXIST **** ') IF(MODE) 11,3,16 C *** WHEN MODE=0 THE COLOR TRANSFORMATION TABLES ARE READ 3 TBND=LOG10(5500.) GBND=2. OPEN(UNIT=12,ERR=49,FILE='uvbylo.data',STATUS='OLD') OPEN(UNIT=13,ERR=49,FILE='uvbyhi.data',STATUS='OLD') READ(12,1000) NT,NG,NFE,NDX READ(12,1001) (TA(I),I=1,NT) READ(12,1002) (G(I),I=1,NG) DO 4 I=1,NT 4 T(I)=LOG10(TA(I)) DO 5 K=1,NFE READ(12,1006) NAMFE,FEH(K),NAMCOL DO 5 J=1,NG READ(12,1003) (BY(I,J,K),I=1,NT) 5 CONTINUE DO 6 K=1,NFE READ(12,1006) NAMFE,FEH(K),NAMCOL DO 6 J=1,NG READ(12,1003) (M1(I,J,K),I=1,NT) 6 CONTINUE DO 7 K=1,NFE READ(12,1006) NAMFE,FEH(K),NAMCOL DO 7 J=1,NG READ(12,1003) (C1(I,J,K),I=1,NT) 7 CONTINUE DO 8 K=1,NFE READ(12,1006) NAMFE,FEH(K),NAMCOL DO 8 J=1,NG READ(12,1003) (BC(I,J,K),I=1,NT) 8 CONTINUE READ(13,1000) NTT,NGG,NFE,NDXX READ(13,1001) (TB(I),I=1,NTT) READ(13,1002) (GG(I),I=1,NGG) DO 9 I=1,NTT 9 TT(I)=LOG10(TB(I)) DO 10 L=1,NDXX DO 10 K=1,NFE READ(13,1006) NAMFE,FEH(K),NAMCOL DO 10 J=1,NGG READ(13,1003) (A(I,J,K,L),I=1,NTT) 10 CONTINUE CLOSE(UNIT=12,STATUS='KEEP') CLOSE(UNIT=13,STATUS='KEEP') C *** WHEN MODE=0 OR MODE=-1, COLOR TRANSFORMATION TABLES ARE C *** CREATED FOR THE INPUT [Fe/H] VALUE USING LINEAR INTERPOLATION 11 DO 12 M=2,NFE K=M-1 IF(FE.LE.FEH(M)) GO TO 13 12 CONTINUE 13 M=K+1 SLOPE=(FE-FEH(K))/(FEH(M)-FEH(K)) DO 14 J=1,NG DO 14 I=1,NT C(I,J)=BY(I,J,K)+SLOPE*(BY(I,J,M)-BY(I,J,K)) D(I,J)=M1(I,J,K)+SLOPE*(M1(I,J,M)-M1(I,J,K)) E(I,J)=C1(I,J,K)+SLOPE*(C1(I,J,M)-C1(I,J,K)) F(I,J)=BC(I,J,K)+SLOPE*(BC(I,J,M)-BC(I,J,K)) 14 CONTINUE DO 15 L=1,NDXX DO 15 J=1,NGG DO 15 I=1,NTT B(I,J,L)=A(I,J,K,L)+SLOPE*(A(I,J,M,L)-A(I,J,K,L)) 15 CONTINUE WRITE(6,1004) FE C *** WHEN MODE=0 OR MODE=-1, CONTROL RETURNS TO THE CALLING C *** PROGRAM ONCE THE [Fe/H] INTERPOLATION IS CARRIED OUT (I.E., NO C *** INTERPOLATIONS ARE PERFORMED FOR INPUT VALUES OF GV AND TEFF) GO TO 50 C *** WHEN MODE=1, INTERPOLATIONS ARE CARRIED OUT FOR THE INPUT VALUES C *** OF [Fe/H], log g, AND log Teff (= FE, GV, AND TEFF IN THE ARGUMENT C *** LIST) 16 IF(TEFF.LE.TBND) GO TO 18 IF(GV.GE.GBND) GO TO 40 IF(TEFF.LE.T(NT)) GO TO 18 C *** EXECUTION HALTS WITH A "STOP30" CODE IF THE INPUT TEMPERATURE IS C *** OUTSIDE THE RANGE OF THE TABLES ON THE HIGH SIDE WRITE(8,1005) GV,TEFF STOP30 C *** THE NEXT SECTION ASSUMES THAT THE LOW-TEMPERATURE TABLES ARE THE C *** RELEVANT ONES TO USE 18 NGM=NG-1 DO 19 I=3,NGM MG=I-2 IF(GV.LE.G(I)) GO TO 20 19 CONTINUE IF(GV.LE.G(NG)) GO TO 20 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log g CONSIDERED IN THE C *** TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT TWO C *** LINES SHOULD BE ACTIVATED C WRITE(8,1005) GV,TEFF C STOP31 20 R(1)=G(MG) R(2)=G(MG+1) R(3)=G(MG+2) R(4)=G(MG+3) CALL LGRAN4(R,AG,GV) NTM=NT-1 DO 25 I=3,NTM MT=I-2 IF(TEFF.LE.T(I)) GO TO 30 25 CONTINUE IF(TEFF.LE.T(NT)) GO TO 30 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log Teff CONSIDERED IN C *** THE TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT C *** TWO LINES SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP32 30 R(1)=T(MT) R(2)=T(MT+1) R(3)=T(MT+2) R(4)=T(MT+3) CALL LGRAN4(R,AT,TEFF) L=MG-1 DO 35 I=1,4 L=L+1 O(I)=AT(1)*C(MT,L)+AT(2)*C(MT+1,L)+AT(3)*C(MT+2,L)+AT(4)*C(MT+3,L) P(I)=AT(1)*D(MT,L)+AT(2)*D(MT+1,L)+AT(3)*D(MT+2,L)+AT(4)*D(MT+3,L) Q(I)=AT(1)*E(MT,L)+AT(2)*E(MT+1,L)+AT(3)*E(MT+2,L)+AT(4)*E(MT+3,L) 35 R(I)=AT(1)*F(MT,L)+AT(2)*F(MT+1,L)+AT(3)*F(MT+2,L)+AT(4)*F(MT+3,L) C *** 4-POINT LAGRANGIAN INTERPOLATION IS USED TO DERIVE THE VALUES OF C *** b-y, m1, c1, AND BC_V CORRESPONDING TO THE INPUT VALUES OF C *** [Fe/H], log g, and log Teff. SBY=AG(1)*O(1)+AG(2)*O(2)+AG(3)*O(3)+AG(4)*O(4) SM1=AG(1)*P(1)+AG(2)*P(2)+AG(3)*P(3)+AG(4)*P(4) SC1=AG(1)*Q(1)+AG(2)*Q(2)+AG(3)*Q(3)+AG(4)*Q(4) BCV=AG(1)*R(1)+AG(2)*R(2)+AG(3)*R(3)+AG(4)*R(4) C *** SPLINE INTERPOLATION IS USED TO FIND THE VALUE OF BC_V AT LOW C *** TEMPERATURES. DO 38 I=1,NT X(I)=T(I) 38 Y(I)=AG(1)*F(I,MG)+AG(2)*F(I,MG+1)+AG(3)*F(I,MG+2)+AG(4)*F(I,MG+3) CALL INTEP(TEFF,BCV,X,Y,NT,IER) GO TO 50 C *** THE NEXT SECTION ASSUMES THAT THE HIGH-TEMPERATURE TABLES ARE THE C *** RELEVANT ONES TO USE. 40 IF(TEFF.LE.TT(NTT)) GO TO 42 WRITE(8,1005) GV,TEFF STOP33 42 NGM=NGG-1 DO 43 I=3,NGM MG=I-2 IF(GV.LE.GG(I)) GO TO 44 43 CONTINUE IF(GV.LE.GG(NGG)) GO TO 44 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log g CONSIDERED IN THE C *** TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT TWO C *** LINES SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP34 44 R(1)=GG(MG) R(2)=GG(MG+1) R(3)=GG(MG+2) R(4)=GG(MG+3) CALL LGRAN4(R,AG,GV) NTM=NTT-1 DO 45 I=3,NTM MT=I-2 IF(TEFF.LE.TT(I)) GO TO 46 45 CONTINUE IF(TEFF.LE.TT(NTT)) GO TO 46 C *** SOME EXTRAPOLATION OUTSIDE THE RANGE OF log Teff CONSIDERED IN C *** THE TABLES IS PERMITTED. TO AVOID ANY EXTRAPOLATIONS, THE NEXT C *** TWO LINES SHOULD SHOULD BE ACTIVATED. C WRITE(8,1005) GV,TEFF C STOP35 46 R(1)=TT(MT) R(2)=TT(MT+1) R(3)=TT(MT+2) R(4)=TT(MT+3) CALL LGRAN4(R,AT,TEFF) L=MG-1 DO 48 I=1,4 L=L+1 DO 48 K=1,NDXX H(I,K)=AT(1)*B(MT,L,K)+AT(2)*B(MT+1,L,K)+AT(3)*B(MT+2,L,K)+ 1 AT(4)*B(MT+3,L,K) 48 CONTINUE C *** 4-POINT LAGRANGIAN INTERPOLATION IS USED TO DERIVE THE VALUES OF C *** b-y, m1, c1, AND BC_V CORRESPONDING TO THE INPUT VALUES OF C *** [Fe/H], log g, and log Teff. SBY=AG(1)*H(1,1)+AG(2)*H(2,1)+AG(3)*H(3,1)+AG(4)*H(4,1) SM1=AG(1)*H(1,2)+AG(2)*H(2,2)+AG(3)*H(3,2)+AG(4)*H(4,2) SC1=AG(1)*H(1,3)+AG(2)*H(2,3)+AG(3)*H(3,3)+AG(4)*H(4,3) BCV=AG(1)*H(1,4)+AG(2)*H(2,4)+AG(3)*H(3,4)+AG(4)*H(4,4) GO TO 50 49 WRITE(6,1007) STOP 50 RETURN END SUBROUTINE LGRAN4(X,A,XX) IMPLICIT DOUBLE PRECISION (A-H, O-Z) SAVE DIMENSION X(4),A(4) R1=(X(1)-X(2))*(X(1)-X(3))*(X(1)-X(4)) R2=(X(2)-X(1))*(X(2)-X(3))*(X(2)-X(4)) R3=(X(3)-X(1))*(X(3)-X(2))*(X(3)-X(4)) R4=(X(4)-X(1))*(X(4)-X(2))*(X(4)-X(3)) A(1)=((XX-X(2))*(XX-X(3))*(XX-X(4)))/R1 A(2)=((XX-X(1))*(XX-X(3))*(XX-X(4)))/R2 A(3)=((XX-X(1))*(XX-X(2))*(XX-X(4)))/R3 A(4)=((XX-X(1))*(XX-X(2))*(XX-X(3)))/R4 RETURN END SUBROUTINE INTEP(XP,P,X,F,N,IER) C *** PURPOSE: To interpolate a function value P for a given argument C XP from a table of N values (X,F). This is a spline interpolation C scheme based on Hermite polynomials. The source is U.S. Airforce C Surveys in Geophysics No 272. C *** USAGE: For random values of XP C CALL INTEP(XP,P,X,F,N,IER) C or after the first call to INTEP with monotonically increasing or C decreasing values of XP consistent with the X vector C CALL EINTEP(XP,P,X,F,N,IER) C DESCRIPTION OF PARAMETERS: C XP - the chosen argument value C P - the resultant interpolated value C X - the vector of independent values C F - the vector of dependent values C N - the number of points in the (X,F) vectors C IER - the resultant error parameter (set = 2 if XP is beyond C either extreme of X; in which case P is set to the value of C F at that extremum) IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 LP1,LP2,L1,L2 DIMENSION F(*),X(*) IER=1 IO=1 IUP=0 IF(X(2).LT.X(1)) IUP=1 N1=N-1 IF((XP.GE.X(N).AND.IUP.EQ.0).OR.(XP.LE.X(N).AND.IUP.EQ.1)) THEN * 5 P=F(N) P=F(N) * GO TO 6 IER=2 RETURN ELSE IF((XP.LE.X(1).AND.IUP.EQ.0).OR. 1 (XP.GE.X(1).AND.IUP.EQ.1)) THEN P=F(1) * 6 IER=2 IER=2 RETURN ENDIF ENTRY EINTEP(XP,P,X,F,N,IER) 8 DO 1 I=IO,N IF(XP.LT.X(I).AND.IUP.EQ.0) GO TO 2 IF(XP.GT.X(I).AND.IUP.EQ.1) GO TO 2 1 CONTINUE P=F(N) IER=2 RETURN * GO TO 5 2 I=I-1 IF(I.EQ.IO-1) GO TO 4 IO=I+1 LP1=1./(X(I)-X(I+1)) LP2=1./(X(I+1)-X(I)) IF(I.EQ.1) FP1=(F(2)-F(1))/(X(2)-X(1)) IF(I.EQ.1) GO TO 3 FP1=(F(I+1)-F(I-1))/(X(I+1)-X(I-1)) 3 IF(I.GE.N1) FP2=(F(N)-F(N-1))/(X(N)-X(N-1)) IF(I.EQ.N1) GO TO 4 FP2=(F(I+2)-F(I))/(X(I+2)-X(I)) 4 XPI1=XP-X(I+1) XPI=XP-X(I) L1=XPI1*LP1 L2=XPI*LP2 P=F(I)*(1.-2.*LP1*XPI)*L1*L1+F(I+1)*(1.-2.*LP2*XPI1)*L2*L2+ 1 FP2*XPI1*L2*L2+FP1*XPI*L1*L1 RETURN END