program mash3iso C This is a version without error propagation or error messages C concerning possible wild extrapolations. (But it runs twice as C fast as the regular mash3 program.) integer nind,nvk,nfe,ng parameter (nind=9,nvk=86,nfe=7,ng=13) C the 8th nind is the temperature. C There are actually only 75 vk points, but don't change the nvk param. C the 9th nind is BC_v real a(nvk,nfe,ng,nind) real grav,feh,teff,theta,clrs(nind) C out-of-bounds flag. Nominal is zero. integer iflag C load data table call readtable(a) C inquire as to input parameters print*, 'enter log g ' read*, grav print*, 'enter [Fe/H]' read*, feh print*, 'enter Teff ' read*, teff print*, ' ' C switch Teff to THETA theta = 5040./teff C get colors for these parameters call teffinterp(a,grav,feh,theta,clrs,iflag) C iflag, which is ignored here, warns of extrapolations. C Nominally zero, it becomes -N if too hot or N if too cool, C where N is the number of point spacings extrapolated. C Obviously, abs(iflag) > 1 is of concern . . . write(*,'(a6,f7.3,a5,f5.3)') 'U-B = ',clrs(1) write(*,'(a6,f7.3,a5,f5.3)') 'B-V = ',clrs(2) write(*,'(a6,f7.3,a5,f5.3)') 'V-R = ',clrs(3) write(*,'(a6,f7.3,a5,f5.3)') 'V-I = ',clrs(4) write(*,'(a6,f7.3,a5,f5.3)') 'J-K = ',clrs(5) write(*,'(a6,f7.3,a5,f5.3)') 'H-K = ',clrs(6) write(*,'(a6,f7.3,a5,f5.3)') 'V-K = ',clrs(7) write(*,'(a6,f7.3,a5,f5.3)') 'BCv = ',clrs(9) stop end C-------------------------------------------------------------------------- subroutine readtable(a) integer nind,nvk,nfe,ng,i parameter (nind=9,nvk=86,nfe=7,ng=13) C the 8th nind is the temperature C the 9th nind is BC_v real a(nvk,nfe,ng,nind),g(ng),fe(nfe) data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 / data fe / -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 / open(unit=19,file='mash2.out',status='old') do ife = 1,nfe do ig = 1,ng do ivk=1,75 read(19,*) x1,x2,xteff,x3, z (a(ivk,ife,ig,k),k=1,7),a(ivk,ife,ig,9) if ( abs(x1 - fe(ife)).gt. 0.01 ) print*, 'bad fe' if ( abs(x2 - g(ig)).gt. 0.01 ) print*, 'bad g' C transform to THETA rather than Teff . . . a(ivk,ife,ig,8) = 5040./xteff end do end do end do return end C---------------------------------------------------------------------------- C subroutine teffinterp subroutine teffinterp(a,grav,feh,theta,clrs,iflag) C given a THETA (=5040/Teff), return colors integer nind,nvk,nfe,ng,iflag parameter (nind=9,nvk=86,nfe=7,ng=13) real grav,feh,theta,clrs(nind) real a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk) C the 8th nind is the temperature, 9th BC_v C local variables integer jg,jfe real c(nvk,nind) data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 / data fe / -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 / C clrs are 1 U-B 2 B-V 3 V-R 4 V-I 5 J-K 6 H-K 7 V-K (8 Teff) C a() is organised: a(nvk,nfe,ng,nind) (9 BCv) C find jfe and jg interpolation corners call locate(g,ng,grav,jg) if (jg.eq.0) jg = 1 if (jg.eq.ng) jg = ng-1 call locate(fe,nfe,feh,jfe) if (jfe.eq.0) jfe=1 if (jfe.eq.nfe) jfe=nfe-1 C fill c array with bilinear-interp results ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe) ) gg = ( grav - g(jg) )/( g(jg+1) - g(jg) ) do ivk=1,75 do ind=1,nind c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind) z + ffe *(1.0-gg)*a(ivk,jfe+1,jg,ind) z + ffe *gg *a(ivk,jfe+1,jg+1,ind) z + (1.0-ffe)*gg *a(ivk,jfe,jg+1,ind) end do end do C find temperature (it's really THETA) and interpolate colors call locate(c(1,8),75,theta,jt) jt0 = jt if (jt.eq.0) jt = 1 if (jt.gt.70) jt = 71 do i=1,7 call polint(c(jt,8),c(jt,i),5,theta,clrs(i),dy) end do call polint(c(jt,8),c(jt,9),5,theta,clrs(9),dy) C if requested temperature is out-of-bounds, use linear interpolation to C extrapolate. Return iflag = int(number of segments beyond the tabulated) if ( jt0 .eq. 0 ) then frac = (theta - c(1,8))/(c(2,8)-c(1,8)) iflag = -1 + int(frac) do i=1,7 clrs(i) = (1.0-frac)*c(1,i) + frac*c(2,i) end do clrs(9) = (1.0-frac)*c(1,9) + frac*c(2,9) end if if ( jt0 .eq. 75 ) then frac = (theta - c(74,8))/(c(75,8)-c(74,8)) iflag = 1 + int(frac) do i=1,7 clrs(i) = (1.0-frac)*c(74,i) + frac*c(75,i) end do clrs(9) = (1.0-frac)*c(74,9) + frac*c(75,9) end if return end C----------------------------------------------------------------------- C subroutine vkinterp subroutine vkinterp(a,grav,feh,theta,clrs) C given a V-K (in clrs(7)), return colors and THETA = 5040/Teff integer nind,nvk,nfe,ng parameter (nind=9,nvk=86,nfe=7,ng=13) real grav,feh,theta,clrs(nind) real a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk) C the 8th nind is the temperature, 9th BC_v C local variables integer jg,jfe real c(nvk,nind) data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 / data fe / -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 / C clrs are 1 U-B 2 B-V 3 V-R 4 V-I 5 J-K 6 H-K 7 V-K (8 Teff) C a() is organised: a(nvk,nfe,ng,nind) (9 BCv) C find jfe and jg interpolation corners call locate(g,ng,grav,jg) if (jg.eq.0) jg = 1 if (jg.eq.ng) jg = ng-1 call locate(fe,nfe,feh,jfe) if (jfe.eq.0) jfe=1 if (jfe.eq.nfe) jfe=nfe-1 C fill c array with bilinear-interp results ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe) ) gg = ( grav - g(jg) )/( g(jg+1) - g(jg) ) do ivk=1,75 do ind=1,nind c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind) z + ffe *(1.0-gg)*a(ivk,jfe+1,jg,ind) z + ffe *gg *a(ivk,jfe+1,jg+1,ind) z + (1.0-ffe)*gg *a(ivk,jfe,jg+1,ind) end do end do C find V-K and interpolate THETA and colors call locate(c(1,7),75,clrs(7),jt) if (jt.eq.0) jt = 1 if (jt.gt.70) jt = 71 do i=1,6 call polint(c(jt,7),c(jt,i),5,clrs(7),clrs(i),dy) end do call polint(c(jt,7),c(jt,8),5,clrs(7),clrs(8),dy) call polint(c(jt,7),c(jt,9),5,clrs(7),clrs(9),dy) theta = clrs(8) return end C----------------------------------------------------------------------- C subroutine viinterp subroutine viinterp(a,grav,feh,theta,clrs) C given a V-I (in clrs(4)), return colors and THETA = 5040/Teff integer nind,nvk,nfe,ng parameter (nind=9,nvk=86,nfe=7,ng=13) real grav,feh,theta,clrs(nind) real a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk) C the 8th nind is the temperature, 9th BCv C local variables integer jg,jfe real c(nvk,nind) data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 / data fe / -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 / C clrs are 1 U-B 2 B-V 3 V-R 4 V-I 5 J-K 6 H-K 7 V-K (8 Teff) C a() is organised: a(nvk,nfe,ng,nind) C find jfe and jg interpolation corners call locate(g,ng,grav,jg) if (jg.eq.0) jg = 1 if (jg.eq.ng) jg = ng-1 call locate(fe,nfe,feh,jfe) if (jfe.eq.0) jfe=1 if (jfe.eq.nfe) jfe=nfe-1 C fill c array with bilinear-interp results ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe) ) gg = ( grav - g(jg) )/( g(jg+1) - g(jg) ) do ivk=1,75 do ind=1,nind c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind) z + ffe *(1.0-gg)*a(ivk,jfe+1,jg,ind) z + ffe *gg *a(ivk,jfe+1,jg+1,ind) z + (1.0-ffe)*gg *a(ivk,jfe,jg+1,ind) end do end do C find V-I and interpolate THETA and colors call locate(c(1,4),75,clrs(4),jt) if (jt.eq.0) jt = 1 if (jt.gt.70) jt = 71 do i=1,3 call polint(c(jt,4),c(jt,i),5,clrs(4),clrs(i),dy) end do do i=5,9 call polint(c(jt,4),c(jt,i),5,clrs(4),clrs(i),dy) end do theta = clrs(8) return end C----------------------------------------------------------------------- C------------------------------------------------------------------------ C subroutine linear C quick linear interpolation subroutine linear(x,y,nxy,xin,yout,iOK) integer nxy,iOK real x(nxy),y(nxy),xin,yout C x and y are nxy long C xin is the nontabulated input x value C yout is the interpolated y guess. C iOK is -3 if X is not in ascending order, -2 if input x is C out-of-bounds by more than 1 xpoint spacing, -1 if out-of-bounds C by less than 1 xpoint spacing, 0 if all is OK C local variables integer j,jl,ju,jm real frac iOK = 0 if ( x(2).lt.x(1) ) then print*, z 'Error. Sub LINEAR. Input X array must be in ascending order.' yout = 0.0 iOK = -3 return end if C locate correct array element by bisection jl=0 ju=nxy+1 10 if (ju-jl.gt.1) then jm = (ju+jl)/2 if ((x(nxy).gt.x(1)).eqv.(xin.gt.x(jm))) then jl=jm else ju=jm end if goto 10 end if j = jl C j is 0 or nxy if xin is off the grid C if off-grid, reset j and set output flag iOK if ( j.eq.0) then if ( xin.lt.(x(1)-(x(2)-x(1))) ) then iOK = -2 else iOK = -1 end if j=1 endif if ( j.eq.nxy) then if ( xin.gt.(x(nxy)+(x(nxy)-x(nxy-1))) ) then iOK = -2 else iOK = -1 end if j = nxy-1 end if C now interpolate/extrapolate frac = (xin - x(j))/(x(j+1)-x(j)) yout = (1.0-frac)*y(j) + frac*y(j+1) return end C end subroutine linear --------------------------------------------- C -----NUMERICAL RECIPES routines: locate and polint SUBROUTINE LOCATE(XX,N,X,J) DIMENSION XX(N) JL=0 JU=N+1 10 IF(JU-JL.GT.1)THEN JM=(JU+JL)/2 IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM)))THEN JL=JM ELSE JU=JM ENDIF GO TO 10 ENDIF J=JL RETURN END SUBROUTINE POLINT(XA,YA,N,X,Y,DY) PARAMETER (NMAX=10) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END