C------------------------------------------------------------------------------ C C Bahcall-Soneira Galaxy model for star counts October 1978 C Visual Band counts with B-V colors C C Provision for transformation to February 1983 C other bands and colors C R.M. Soneira c c Includes reddening as well as obscuration. R. M. Soneira/ c J. N. Bahcall c (about 6/85) c c c Includes cutoff in disk giants at Mabs = -1.5 (B-V = 1.7) Bahcall c (5/86) c Uses IMSL subroutines for spline fitting of color-magnitude c diagrams. c C------------------------------------------------------------------------------ c C The model is based on the following papers: C 1. Ap. J. Suppl. 44, 73 (1980) - Original model paper C 2. Ap. J. Letters 238, L17 (1980) - Principal results from above C 3. Ap. J. 246, 122 (1981) - Bright stars and disk giants C 4. Ap. J. Suppl. 47, 401 (1981) - Tabulated counts; halo counts C 5. Ap. J. 255, 181 (1982) - Luminosity function in R,I,J,K bands C 6. Ap. J. 265, 730 (1983) - Galactic spheroid in detail C 7. Ap. J. 272, 627 (1983) - Giants in the spheroid C 8. Ap. J.,Suppl. 55, 67 (1984) - Comparisons with 5 Fields C 9. IAU Colloquium 76, p. 259 (1983)-Discussion and comparison c 10. Ap. J. 299, December 15(1984) -The Basel Catalog c 11 Ann.Rev.Astr.Astrp.24,577 (1986)-Summary/Next Steps/Open Questions c C------------------------------------------------------------------------------ c c c *********************************************************************** c c For the purpose of this program only, IMSL has authorized the reprint c of selected source code taken from the IMSL library. This source code c is subject to copyright protection and is considered to be confidential c and proprietary information of IMSL. Use, duplication, or distribution c of the IMSL source code in any form is strictly prophibited. c c c ************************************************************************ C C The parameter values specified are for C The Bahcall-Soneira Standard Galaxy Model for the Visual band. C C Those parameters which change with band are marked with the symbol ### C See paper 5 for parameter values for the R, I, J, K bands. C See papers 4 and 5 for transformation methods to other bands. C C It is possible to keep the absolute magnitude distribution in space C in the Visual band and transform only the apparent magnitudes C to any band which can be expressed in terms of V and B-V. C Colors can also be transformed from B-V into any other system which C can be expressed in terms of B-V and/or absolute Visual magnitude. C The place to do this is in the subroutine BVC which C knows both V and B-V. C C Note: If you are transforming only the apparent magnitudes, C then use parameters for V absolute magnitudes throughout. C------------------------------------------------------------------------------ C The model has been compared to data down to 20 deg galactic latitude. C This program has been tested to only 10 deg galactic latitude. C Please exercise caution below these limits. C------------------------------------------------------------------------------ C Input/Output: C unit=5 for terminal read for program queries C unit=6 for terminal write for program replies to queries C unit=8 for terminal/disk file output table of counts and colors C unit=10 output table to be used for plotting color distribution C (column output format is: C color total disk spheroid (in counts per color bin) C see Subroutine COLOR) C------------------------------------------------------------------------------ C Fortran usage complications: C This program is written primarily in a simple subset of Fortran'66 C in order to minimize problems in portability. (A few extensions C to Fortran'66 are used that should be available on most C compilers: IMPLICIT statement, dimensioning in a type statement, C literals inside apostrophes, END and ERR branches inside a C READ statement. C There are two exceptions which are likely to cause some users problems: C Many variables in this program have names that are longer C than 6 characters for clarity and auto documentation. C Some Fortran compilers ignore characters beyond the sixth C in variable names and produce a warning message for each occurrence. C All variable names in this program are unique in the first C 6 characters and therefore the warnings can be ignored. C All STOP statements in this program have a comment field that C prints a message at the terminal to indicate the origin of the C stop (normal or abnormal). Some Fortran compilers do not C allow comments following a stop, but do allow stop number ids. C If this is a problem, change all the comments to unique C numbers, for example: C STOP '*** color convolution too large' C is replaced by C STOP 3 C------------------------------------------------------------------------------ C This program has been tested on: C an IBM/370 CMS system with the Fortran'66 G, H, and VS compilers, C a Digital VAX/VMS system running the DEC Fortran'77 compiler, C and on the Berkeley Unix version 4.1 & 4.2 f77 Fortran'77 compilers. C C note for f77 Unix users: some versions of the f77 compiler have C a bug which treats variables longer than 6 characters as C an error condition and not a warning condition. This C causes compilation to be terminated. In order to get C around this problem use the program version model6.f on C this distribution tape, or if you have obtained this program C via other channels, use the error listing produced by the C compiler and an editor to truncate the offending variables. C As mentioned above, the variables are unique in the first C six characters. C C Another alternative is to fix the bug as follows: C change the second of the lines below from /usr/src/cmd/f77/lex.c C sprintf(buff, "name %s too long, truncated to %d",token, VL); C err(buff); C to C sprintf(buff, "name %s too long, truncated to %d",token, VL); C warn(buff); C then run "make" in /usr/src/cmd/f77 to recompile the compiler C------------------------------------------------------------------------------ C IMPLICIT REAL (A-H,L,M,O-Z) C dimension arrays for number-magnitude counts REAL XNM(100), XNMS(100) REAL LFD(1000), LFS(1000), NS DOUBLE PRECISION BLANK C spheroid color-magnitude diagram ids DOUBLE PRECISION GIANTBR(7) C COMMON /LF/ NS,MS,A,B,D,ABD,MDIM,MBR,MTURN,BLUESHIFT,MSUN COMMON /PARM/ R0,A0,ABSM,DM,DMA,MLIM,ML1,ML2,ASIZE,L2,B2,DRDISK,S DATA GIANTBR/'M67 ', 'M3 ', 'M5 ', 'M13 ', 1 'M15 ', 'M92 ', '47Tuc'/ DATA XNM/100*0.0/, XNMS/100*0.0/ C C------------------------------------------------------------------------------ C This program will produce a large number of floating point C exponent underflows, which are in general harmless. C It is best to run with underflows masked off. C On VMS systems this is often the system default, if not it C can be set by calling the VMS "Condition Handling Facility" C C ON IBM systems the following two subroutine calls will accomplish this. C suppress underflow and tangent overflow messages <--- for IBM only C CALL ERRSET (208,256,-1,1) C CALL ERRSET (259,256,-1,1) C------------------------------------------------------------------------------ c c Use these Open statements only on a VAX c Open(8,file = 'redden.out', status = 'unknown') Open(10,file = 'redden.colortable', status = 'unknown') Write(10,137) 137 Format(/,1x,'Columns are: Color, Total, Disk, Spheroid.', $/,2x, 'Last three are in counts per color bin.',/, $2x, 'The counts are in number per square degree.',//) Write(10,138) 138 Format(5x,'Color Total Disk Spheroid',/) c c The color table normalization, BVSC, was checked on 7/29/91. C PIE = 3.141593 C = PIE/180.0 C viewing size (square degrees) ASIZE = 1.0 C convert to steradians OMEGA = ASIZE * C*C C C------------------------------------------------------------------------------ C Default Parameters C C dimmest apparent magnitude for calculation MLIM = 31 C apparent magnitude bin size DM = 0.5 C luminosity function dim-end cutoff ### MDIM = +16.5 C luminosity function bright-end cutoff for disk ### MBR = -6 C luminosity function bright-end cutoff for spheroid (CUTOFF > MBR) ### CUTOFF = -3 C luminosity function parameters. See Papers 1, 5, 6, 7. C Values are for the Visual band. See paper 5 for other bands. ### NS = 4.03E-3 MS = 1.28 MC = 15.0 A = 0.74 B = 0.04 D = 3.40 ABD= (A-B)/D C absolute magnitude integration step size DMA = 0.05 C default color interval for printout DBVPR = 0.10 C default color bin size for scaling color plot (magnitudes) BVPLBIN = 0.20 C default color interval for plotting DBVPL = 0.01 C round to submultiple of print interval ICOMPR = DBVPR/DBVPL + 0.5 IF (ICOMPR .LT. 1) ICOMPR = 1 DBVPL = DBVPR/ICOMPR BVSC = BVPLBIN/DBVPL c c Set BVSC = 1.0 to get numers per square degree. Use BVSC = 1./3600 c to get normalization per square arc minute. 9/29/91. BVSC = 1.00 c C default color error in magnitudes SBV = 0.10 C internal color interval for calculation before convolution DBV = SBV/20 C round to submultiple of print interval ICOMPL = DBVPL/DBV + 0.5 IF (ICOMPL .LT. 1) ICOMPL = 1 DBV = DBVPL/ICOMPL C default maximum B-V color BVMAX = 4 c BVMAX = 3 c ratio of obscuration to reddening Av/E(B-V) = 3 AEBAND = 3.0 C default giant branch for spheroid ( M13 ). See array GIANTBR at top. IGIANTBR = 4 C disk / spheroid star density in solar neighborhood DENCON = 500.0 C default apparent magnitude interval for color distribution ML1 = 0 ML2 = 20 C default angular direction L2 = 0 B2 = 90 C distance of Sun from Galactic center (pc) R0 = 8000 C scale height for absorption (pc) A0 = 100 C spheroid axis ratio E = Rez/Rex and Rey = Rex E = 0.8 C de Vaucouleurs radius for spheroid (kpc). Rex * Rez = Re*Re = REK*REK REK = 2.67 C default components for calculation [1=disk 2=spheroid 3=both] ICOMP = 3 C default plotting [-1=no +1=yes] IPLOT = +1 C C------------------------------------------------------------------------------ C prompt for input values C C The program will print out the default values for several parameters. C To override a default value, type in a new value in the appropriate C column position. Since the input format is F6.2, to override, e.g., C the third parameter, skip 12 spaces and type the desired value in C the next six spaces, using standard fortran rules for input. C Inputing blanks or zeros in a field results in the printed default C parameter value being adopted. Entering a carriage return is C equivalent to entering blanks for all subsequent parameters C on the line. C C------------------------------------------------------------------------------ 10 CONTINUE 11 WRITE (6,1) L2,B2,IPLOT,ICOMP,E,REK 1 FORMAT (//' Enter in [f6.2] L',T32,'B',T40,'plot?', 1 T49,'components',T65,'e',T72,'Re(kpc)'/ 2 ' defaults are:',F10.2,F10.2,I8,I11,F14.2,F11.2) READ (5,2,END=1000,ERR=11) L2I,B2I,PLOTI,COMPI,EI,REKI 2 FORMAT (6F6.2) C inputing a blank, zero or a carriage return means previous value IF (L2I .NE. 0) L2 = L2I IF (B2I .NE. 0) B2 = B2I IF (EI .NE. 0) E = EI IF (REKI .NE. 0) REK = REKI C convert to kpc RE = REK*1E3 IF (PLOTI .GT. 0) IPLOT = +1 IF (PLOTI .LT. 0) IPLOT = -1 IF (COMPI .GT. 0) ICOMP = COMPI + 0.5 WRITE (6,101) L2,B2,IPLOT,ICOMP,E,REK 101 FORMAT (' Adopted values:',F8.2,F10.2,I8,I11,F14.2,F11.2/) C C Sandage obscuration model for Visual band ### ABSM = 0.165 * ( 1.192-TAN(ABS(B2)*C) ) / SIN(ABS(B2)*C) C zero obscuration for latitudes > 50 deg IF (ABS(B2) .GE. 50) ABSM = 0 c c Insert here the Burstein-Heiles values for ABSM on right-hand c side of the following equation. Can use this instead of the c Sandage model. c c Absm = Burstein-Heiles value [3*E(B-V)] C 12 WRITE (6,3) ML1,ML2,MDIM,MBR,DMA 3 FORMAT (/' Enter in [f6.2]',T24,'m1',T31,'m2',T43,'Mdim', 1 T51,'Mbr',T64,'dM'/ 2 ' defaults are:',3(F12.2,F8.2)) READ (5,4,END=500,ERR=12) ML1I,ML2I,MDIMI,MBRI,DMAI 4 FORMAT (5F6.2) C inputing a blank, zero or a carriage return means previous value IF (ML1I .NE. 0) ML1 = ML1I IF (ML2I .NE. 0) ML2 = ML2I IF (MDIMI .NE. 0) MDIM = MDIMI IF (MBRI .NE. 0) MBR = MBRI IF (DMAI .NE. 0) DMA = DMAI 500 CONTINUE WRITE (6,103) ML1,ML2,MDIM,MBR,DMA 103 FORMAT (' Adopted values:',F10.2,F8.2,F12.2,F8.2,F12.2/) C C enter spheroid giant branch id 600 WRITE (6,7) GIANTBR(IGIANTBR),(I,GIANTBR(I),I=1,7) 7 FORMAT (/' Default Spheroid Giant branch=',A5, 1 5X/' Enter: [ ',7(I1,'=',A5),' ]') READ (5,8,END=1000,ERR=600) IRGIANTBRI 8 FORMAT (I1) C inputing a blank, zero or a carriage return means previous value IF (IRGIANTBRI .LT. 0 .OR. IRGIANTBRI .GT. 7) GO TO 600 IF (IRGIANTBRI .NE. 0) IGIANTBR = IRGIANTBRI WRITE (6,107) GIANTBR(IGIANTBR) 107 FORMAT (' Adopted Spheroid Giant branch=',A5/) C C------------------------------------------------------------------------------ C begin computation C C set up luminosity function CALL LUMFUN (LFD,LFS) C initialize color routine CALL COLOR (1,GIANTBR,IGIANTBR,IPLOT, 1 BVSC,DBVPR,ICOMPL,ICOMPR,SBV,BVMAX,DBV) C C calculate disk component c c add reddening AEBAND. c IF (ICOMP .NE. 1 .AND. ICOMP .NE. 3) GO TO 60 CALL DISK (XNM,LFD,OMEGA,AEBAND) C print out disk counts CALL PRINT (1,XNM,DM,MLIM,XNMS,L2,B2,ABSM) C C calculate spheroid component 60 IF (ICOMP .NE. 2 .AND. ICOMP .NE. 3) GO TO 70 CALL SPHER (CUTOFF,XNM,LFS,OMEGA,E,RE,DENCON,AEBAND) C print out spheroid counts CALL PRINT (2,XNM,DM,MLIM,XNMS,L2,B2,ABSM) C C now do totals 70 WRITE (8,9) MBR,MDIM,L2,B2,ABSM,ASIZE 9 FORMAT ('TOTAL Mbr=',F5.2,1X,'Mdim=',F5.2, 1 2X,'L=',F6.2, 2 1X,'B=',F5.1,1X,'Obsc=',F5.2,' mag',1X, 3 'Area =',F5.2,' sq deg'/) C print out totals CALL PRINT (3,XNM,DM,MLIM,XNMS,L2,B2,ABSM) C print out color distribution CALL COLOR (2,GIANTBR,IGIANTBR,IPLOT, 1 BVSC,DBVPR,ICOMPL,ICOMPR,SBV,BVMAX,DBV) C 1000 CONTINUE C STOP '*** End of Galaxy Model ***' END SUBROUTINE DISK (XNM,LFD,OMEGA,AEBAND) C compute disk component counts and colors IMPLICIT REAL (A-H,L,M,O-Z) REAL XNM(100), NS REAL LFD(1000) REAL F(1000) C COMMON /LF/ NS,MS,A,B,D,ABD,MDIM,MBR,MTURN,BLUESHIFT,MSUN COMMON /PARM/ R0,A0,ABSM,DM,DMA,MLIM,ML1,ML2,ASIZE,L2,B2,DRDISK,S C PIE = 3.141593 SINB = SIN(ABS(B2)*PIE/180.0) COSB = COS(B2*PIE/180.0) COSL = COS(L2*PIE/180.0) COSBL = COSB*COSL COS2B = COSB*COSB C C scale length in plane (pc) PSH = 3500 C giants scale height (pc) GSH = 250 C C convergence factor for terminating integration CFAC = 1.0E-6 IMAX = MLIM/DM NLF = (MDIM-MBR)/DMA + 1 C C specify fraction of stars on main sequence C relation is for the Visual Band. See paper 5 for other bands ### DO 40 I = 1, NLF MABS = MBR + (I-1)*DMA F(I) = 0.44 * EXP( 1.5E-4 * (MABS+8)**3.5 ) IF (F(I) .GT. 1) F(I) = 1 40 CONTINUE C C distance step size in r direction (pc) DR = 25.0 C maximum height for integration (pc) RMAX = 8000/SINB NR = RMAX/DR C minimum height for integration (pc) RMIN = 0.0 C C------------------------------------------------------------------------------ C integrate along (l,b) direction C TOT = 0 DO 70 J = 1, NR R = J*DR IF (R .LT. RMIN) GO TO 70 C calculate total absorption along line of sight for this distance ABSMAG = ABSM*(1 - EXP(-SINB/A0*R) ) c c calculate reddening for this direction. c Reddening = absmag/aeband c c End of reddening correction c Z = R*SINB X = SQRT(R0*R0 + R*R*COS2B - 2*R*R0*COSBL) EXPX = EXP(- (X-R0)/PSH ) VOL = R*R*OMEGA*DR*DMA ZNM = 0 C C step through luminosity function DO 50 I = 1, NLF MABS = MBR + (I-1)*DMA C calculate apparent magnitude M = MABS + 5.0*ALOG10( R/10 ) + ABSMAG IF (M .GT. MLIM) GO TO 60 C C specify scale heights for main sequence stars C Relation is for the Visual Band ### SH = 82.5*MABS - 97.5 IF (SH .LT. 90.0) SH = 90.0 IF (SH .GT. 325.0) SH = 325.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The variation of disk scale height with absolute magnitude C C in bands other than V is obtained by transforming the two C C characteristic absolute magnitudes that define the transition C C region between small and large scale heights (M=2.3 and M=5.1) C C into other bands using the references cited in papers 1,3,4,5. C C For the B,V,R,I,J,K Johnson bands we get: C C SH = 67.1*MABS - 71.0 For B C C SH = 82.5*MABS - 97.5 For V C C SH = 94.5*MABS - 108 For R C C SH = 106.0*MABS - 124 For I C C SH = 115.2*MABS - 135 For J C C SH = 131.3*MABS - 153 For K C C If you are transforming only the apparent magnitude, then use C C the relation for V absolute magnitudes. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C main sequence and giant fractions FMS = F(I) FG = 1 - FMS c c Cut off disk giants at B-V = 1.7 mag, just as in the globular c clusters. This cutoff implemented in Export programs after May, c 1986. c If(Mabs.LT. -1.5) FG = 0.0 C contribution by this volume element to counts DNMG = VOL*EXPX*LFD(I) * FG*EXP(-Z/GSH) DNMMS = VOL*EXPX*LFD(I) * FMS*EXP(-Z/ SH) DNM = DNMG + DNMMS C total contribution for this step ZNM = ZNM + DNM C C look up colors and perform any desired band transformations CALL BVC (1,MABS,DNMG,DNMMS,DNM,M,MG,MMS,Reddening) C C apparent magnitude bin for giants INDEX = MG/DM IF (INDEX .LT. 1) INDEX = 1 IF (INDEX .GT. IMAX) INDEX = IMAX XNM(INDEX) = XNM(INDEX) + DNMG C C apparent magnitude bin for main sequence stars INDEX = MMS/DM IF (INDEX .LT. 1) INDEX = 1 IF (INDEX .GT. IMAX) INDEX = IMAX XNM(INDEX) = XNM(INDEX) + DNMMS C 50 CONTINUE C 60 CONTINUE C C trivial contribution this distance step? if yes, then done. TOT = TOT + ZNM IF (ZNM .LT. CFAC*TOT) GO TO 80 70 CONTINUE 80 CONTINUE C C print header IRMAX = R + 0.5 WRITE (8,1) MBR,MDIM,IRMAX,L2,B2,ABSM 1 FORMAT ('DISK Mbr=',F5.2,2X,'Mdim=',F5.2, 1 2X,'Rmax=',I6,'pc',2X,'L=',F6.2, 2 1X,'B=',F5.1,2X,'Obsc=',F5.2,' mag'/) RETURN END SUBROUTINE SPHER (CUTOFF,XNM,LFS,OMEGA,E,RE,DENCON,Aeband) C compute spheroid component counts and colors C IMPLICIT REAL (A-H,L,M,O-Z) REAL XNM(100), NS REAL LFS(1000) C COMMON /LF/ NS,MS,A,B,D,ABD,MDIM,MBR,MTURN,BLUESHIFT,MSUN COMMON /PARM/ R0,A0,ABSM,DM,DMA,MLIM,ML1,ML2,ASIZE,L2,B2,DRDISK,S C PIE = 3.141593 SINB = SIN(ABS(B2)*PIE/180.0) COSB = COS(B2*PIE/180.0) COSL = COS(L2*PIE/180.0) COSBL = COSB*COSL C C compute spheroid eccentric factor EC = E RCF = SQRT( COSB*COSB + SINB*SINB/EC/EC ) RE = RE / SQRT(E) REX = RE/1.0E3 REZ = RE*E/1.0E3 C C convergence factor for terminating integration CFAC = 1.0E-6 IMAX = MLIM/DM NLF = (MDIM-MBR)/DMA + 1 C de Vaucouleurs density function parameters BD = 7.66925 C1 = 4.6132E-2/2.0 C2 = PIE/(8*BD) C calculate local value of density RJ = (R0/RE)**0.25 DENR0 = C1*EXP(-BD*RJ) / (RJ*RJ*RJ) * SQRT(C2/RJ) C increment for logarithmic distance integration intervals S = 1.05 X = S - 1 C initial distance (pc) R = 1.0 C maximum distance for integration (pc) RMAX = 1E6 NR = ALOG10(RMAX) / ALOG10(S) C C------------------------------------------------------------------------------ C integrate along (l,b) direction C TOT = 0 DO 70 J = 1, NR R = R*S DR = R*X C calculate total absorption along line of sight for this distance c cDHS prevent underflows by checking the range of -sinb/a0*r being given to exp cDHS modified by dhs 4-oct-1995 to only do exp(-x) for x < 87 cDHS 87 was chosen by testing exp until no underflows occurred on a Sun 4. CDHS old code was: ABSMAG = ABSM*(1 - EXP(-SINB/A0*R) ) c DHSTEST=SINB/A0*R IF(DHSTEST .LT. 87 ) THEN ABSMAG = ABSM*(1 - EXP(-SINB/A0*R) ) ELSE ABSMAG = ABSM ENDIF c c calculate reddening for this distance Reddening = absmag/Aeband C C distance of volume element from center of Galaxy REC = R*RCF RC2 = R0*R0 + REC*REC - 2*R*R0*COSBL RC = SQRT(RC2) VOL = R*R*OMEGA*DR*DMA RJ = (RC/RE)**0.25 C r**1/4 law density DEN = C1*EXP(-BD*RJ) / (RJ*RJ*RJ) * SQRT(C2/RJ) C scale density to local known spheroid star density DEN = DEN / DENR0 / DENCON ZNM = 0 C C step through luminosity function DO 50 I = 1, NLF MABS = MBR + (I-1)*DMA IF (MABS .LT. CUTOFF) GO TO 50 C calculate apparent magnitude M = MABS + 5.0*ALOG10( R/10 ) + ABSMAG IF (M .GT. MLIM) GO TO 60 DNM = VOL*DEN*LFS(I) C total contribution for this distance step ZNM = ZNM + DNM C C main sequence and giant fractions FMS = 1 IF (MABS .LT. MTURN) FMS = 0 FG = 1 - FMS DNMG = DNM * FG DNMMS = DNM * FMS C C look up colors and perform any desired band transformations CALL BVC (2,MABS,DNMG,DNMMS,DNM,M,MG,MMS,Reddening) C C apparent magnitude bin for giants INDEX = MG/DM IF (INDEX .LT. 1) INDEX = 1 IF (INDEX .GT. IMAX) INDEX = IMAX XNM(INDEX) = XNM(INDEX) + DNMG C C apparent magnitude bin for main sequence stars INDEX = MMS/DM IF (INDEX .LT. 1) INDEX = 1 IF (INDEX .GT. IMAX) INDEX = IMAX XNM(INDEX) = XNM(INDEX) + DNMMS 50 CONTINUE C 60 CONTINUE C C trivial contribution this distance step? if yes, then done. TOT = TOT + ZNM IF (ZNM .LT. CFAC*TOT) GO TO 80 C 70 CONTINUE 80 CONTINUE C C print header IRMAX = R/1E3 + 0.5 WRITE (8,1) REX,REZ,IRMAX,E,L2,B2,ABSM 1 FORMAT ('SPHEROID Rex=',F4.2,1X,'Rez=',F4.2, 1 1X,'Rmax=',I4,'kpc',1X,'e=',F4.2,2X,'L=',F6.2, 2 1X,'B=',F5.1,1X,'Obsc=',F5.2,' mag'/) RETURN END SUBROUTINE COLOR (K,GIANTBR,IGIANTBR,IPLOT, 1 BVSC,DBVPR,ICOMPL,ICOMPR,SBV,BVMAX,DBVSET) IMPLICIT REAL (A-H,L,M,O-Z) DOUBLE PRECISION GIANTBR(7) REAL G (500) C COMMON /LF/ NS,MS,A,B,D,ABD,MDIM,MBR,MTURN,BLUESHIFT,MSUN COMMON /PARM/ R0,A0,ABSM,DM,DMA,MLIM,ML1,ML2,ASIZE,L2,B2,DRDISK,S COMMON /BVDIS/ BVMS(1000),BVDISKG(1000),BVSPHERG(1000) COMMON /DBVC/ XBV(2000,3), TBV(3), XGB(3), BVB(3), DBV, IBVMAX C C initialize variables IF (K .GT. 1) GO TO 20 PIE = 3.141593 DBVPL = DBVPR/ICOMPR DBV = DBVPL/ICOMPL IBVMAX = (BVMAX + 1)/DBV + 1 IBVPLMAX = (BVMAX + 1)/DBVPL + 1 IF (IBVMAX .GT. 2000) STOP '*** color bin size too small ***' DBV = DBVSET DO 15 J = 1, 3 TBV(J) = 0 BVB(J) = 0 XGB(J) = 0 DO 10 I = 1, 2000 XBV(I,J) = 0 10 CONTINUE 15 CONTINUE C C fill in main sequence and giant color arrays NLF = (MDIM - MBR)/DMA + 1 DO 30 I = 1, NLF M = MBR + (I-1)*DMA C main sequence BVMS(I) = MAINSEQ(M) C disk giant branch BVDISKG(I) = DISKM67(M) 30 CONTINUE C DO 40 I = 1, NLF M = MBR + (I-1)*DMA C spheroid giant branch selection IF (IGIANTBR .EQ. 1) BVSPHERG(I) = M67 (M) IF (IGIANTBR .EQ. 2) BVSPHERG(I) = M3 (M) IF (IGIANTBR .EQ. 3) BVSPHERG(I) = M5 (M) IF (IGIANTBR .EQ. 4) BVSPHERG(I) = M13 (M) IF (IGIANTBR .EQ. 5) BVSPHERG(I) = M15 (M) IF (IGIANTBR .EQ. 6) BVSPHERG(I) = M92 (M) IF (IGIANTBR .EQ. 7) BVSPHERG(I) = TUC47 (M) 40 CONTINUE C C initialize gaussian convolution function C specify color dispersion CD = SBV IF (CD .LT. DBV) CD = DBV C calculate to 4 sigma NBVCMAX = 4*CD/DBV + 1 IF (NBVCMAX .GT. 500) STOP '*** color convolution too large' DO 50 I = 1, NBVCMAX CC = (I-1)*DBV G(I) = 1.0/SQRT(2.0*PIE)/CD*EXP(-CC*CC/(2*CD*CD))*DBV 50 CONTINUE G(1) = G(1)/2 C RETURN C C------------------------------------------------------------------------------ C output section C C convolve color distribution with gaussian color error 20 CALL CONVOL (XBV,IBVMAX,1,G,NBVCMAX,ICOMPL,IBVPLMAX) CALL CONVOL (XBV,IBVMAX,2,G,NBVCMAX,ICOMPL,IBVPLMAX) C C compute various means BVBARD = BVB(1)/(TBV(1) + 1.0E-10) BVBARS = BVB(2)/(TBV(2) + 1.0E-10) BVBAR = BVB(3)/(TBV(3) + 1.0E-10) XGB(1) = XGB(1)/(TBV(1) + 1.0E-10) XGB(2) = XGB(2)/(TBV(2) + 1.0E-10) XGB(3) = XGB(3)/(TBV(3) + 1.0E-10) FD = TBV(1)/(TBV(3) + 1.0E-10) FS = TBV(2)/(TBV(3) + 1.0E-10) C WRITE (8,1) ML1,ML2,SBV 1 FORMAT ('1[B-V] Distribution for m =',F6.2,' -',F6.2, 1 2X,'Color noise =',F6.3,' mag'/) WRITE (8,12) GIANTBR(IGIANTBR),MTURN,BLUESHIFT 12 FORMAT ('0Spheroid giant branch H.R. diagram = ',A5/ 6 ' Spheroid main sequence turn off =',F5.2,' mag'/ 7 ' Spheroid metallicity color excess =',F5.2,' mag'/) C WRITE (8,2) TBV(3),TBV(1),TBV(2) 2 FORMAT (/10X,'Stars',10X,'=',F13.3,3X,2F12.3) WRITE (8,3) BVBAR,BVBARD,BVBARS 3 FORMAT (10X,'Mean colors',4X,'=',F13.3,3X,2F12.3,9X,7F12.3) WRITE (8,4) FD,FS 4 FORMAT (10X,'Star fraction =',13X,3X,2F12.3) WRITE (8,5) XGB(3),XGB(1),XGB(2) 5 FORMAT (10X,'Giant fraction =',F13.3,3X,2F12.3//) C WRITE (8,6) 6 FORMAT (/' C1',T16,'C2',T25,'C',T35,'Total', 1 T51,'Disk',T60,'Spheroid'/) C C color computation and printout FDISKPR = 0 FSPHRPR = 0 DO 35 I = 1, IBVPLMAX BV2 = I*DBVPL - 1 BV1 = BV2 - DBVPL BV12 = (BV1 + BV2)/2 C FDISK = XBV (I,1) FSPHER = XBV (I,2) FREQ = FDISK + FSPHER FS = FSPHER/(FREQ + 1E-30) C color plotting output. Scale frequencies to desired color bin interval FSC = FREQ * BVSC FDSC = FDISK * BVSC FSSC = FSPHER * BVSC IF (IPLOT .EQ. +1 .AND. I .NE. 1 .AND. I .NE. IBVPLMAX) 1 WRITE (10,9) BV12,FSC,FDSC,FSSC 9 FORMAT (F10.4,10E12.4) C C compress to printing interval IPRPL = I - (I-1)/ICOMPR*ICOMPR IF (IPRPL .EQ. 1) BV1PR = BV1 FDISKPR = FDISKPR + FDISK FSPHRPR = FSPHRPR + FSPHER C C skip printout unless binning is complete IF (IPRPL .NE. ICOMPR) GO TO 35 BV2PR = BV2 BV12PR = (BV1PR + BV2PR)/2 FREQPR = FDISKPR + FSPHRPR FSPR = FSPHRPR/(FREQPR + 1E-30) C only print if non-trivial value IF (FREQPR/TBV(3) .LT. 0.0001) GO TO 34 WRITE (8,7) BV1PR,BV2PR,BV12PR,FREQPR,FDISKPR,FSPHRPR 7 FORMAT (F6.2,' to',F7.2,F8.2,F14.4,3X,2F12.4,2F10.3) C C reset print totals to zero 34 FDISKPR = 0 FSPHRPR = 0 35 CONTINUE C 1000 RETURN END SUBROUTINE CONVOL (F,N,INDEX,G,M,ICCOMP,NOUT) C this routine convolves function F with function G and C compresses the the final result REAL F(2000,3), CON(2000), G(M) C C zero convolution output array DO 50 I = 1, NOUT CON(I) = 0 50 CONTINUE C DO 90 I = 1, N DO 80 J = 1, M K1 = I + J - 1 K2 = I - J + 1 IF (K1 .GT. N) GO TO 75 CON( (K1-1)/ICCOMP + 1) = CON((K1-1)/ICCOMP + 1) + 1 G(J) * F(I,INDEX) 75 IF (K2 .LT. 1) GO TO 80 CON( (K2-1)/ICCOMP + 1) = CON((K2-1)/ICCOMP + 1) + 1 G(J) * F(I,INDEX) 80 CONTINUE 90 CONTINUE C C copy convolved function back to F DO 100 I = 1, NOUT F(I,INDEX) = CON(I) 100 CONTINUE C RETURN END SUBROUTINE BVC (ID,M,FGI,FMSI,FI,MAP,MAPG,MAPMS,reddening) C convert absolute magnitudes to colors C and perform band and color transformations IMPLICIT REAL (A-H,L,M,O-Z) REAL M COMMON /LF/ NS,MS,A,B,D,ABD,MDIM,MBR,MTURN,BLUESHIFT,MSUN COMMON /PARM/ R0,A0,ABSM,DM,DMA,MLIM,ML1,ML2,ASIZE,L2,B2,DRDISK,S COMMON /BVDIS/ BVMS(1000),BVDISKG(1000),BVSPHERG(1000) COMMON /DBVC/ XBV(2000,3), TBV(3), XGB(3), BVB(3), DBV, IBVMAX C C get absolute magnitude index INDEX = (M-MBR)/DMA + 1 C look up main sequence colors BVMSI = BVMS(INDEX) C IF (ID .NE. 1) GO TO 10 C look up disk giant colors BVGI = BVDISKG(INDEX) GO TO 20 C correct for subdwarf branch 10 BVMSI = BVMSI - BLUESHIFT C look up spheroid giant colors BVGI = BVSPHERG(INDEX) C C color band conversion 20 MAPG = MAP + 0.0 * BVGI MAPMS = MAP + 0.0 * BVMSI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C to convert to an apparent magnitude band other than V C C change the above coefficients from zero to the needed values. C C All input and output apparent magnitudes then refer to this new band. C C In general both coefficients should be equal. C C The relation used is pure linear. However, it can be made piecewise C C linear, quadratic, or a nonlinear function. C C C C MAP is the Visual band apparent magnitude C C BVGI is the B-V color for giants C C BVMSI is the B-V color for main sequence stars C C MAPG is the apparent magnitude of giants in the new band C C MAPMS is the apparent magnitude of main sequence stars in the new band C C C C Example: C C B band: m(B) = m(V) + [B-V] C C MAPG = MAP + 1.0 * BVGI C C MAPMS = MAP + 1.0 * BVMSI C C J band: m(J) = m(V) + 0.77[B-V] C C MAPG = MAP + 0.77 * BVGI C C MAPMS = MAP + 0.77 * BVMSI C c g band: m(g) = m(V) + 0.41[B-V] -0.19 C c c The relation for the g-band Gunn filter is given by Kent (1985), c PASP, 97, pg. 165. c c New data for subdwarf color relations is given in the paper by c S. G. Ryan, A. J. 98, 1693 (1989). c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C check if apparent magnitudes are outside the requested range FG = FGI IF (MAPG .LT. ML1 .OR. MAPG .GT. ML2) FG = 0 FMS = FMSI IF (MAPMS .LT. ML1 .OR. MAPMS .GT. ML2) FMS = 0 C C color conversion BVMSI = 1.00 * BVMSI BVGI = 1.00 * BVGI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C to convert to a color other than B-V, C C change the above coefficients from one to the needed values. C C In general both coefficients should be equal. C C The relation used is pure linear. However, it can be made piecewise C C linear, quadratic, or a nonlinear function. C C All input and output colors then refer to this new band. C C C C Example: C C J-F C C BVMSI = 1.17 * BVMSI C C BVGI = 1.17 * BVGI C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c c Add reddening c BVMSI = BVMSI + REDDENING BVGI = BVGI + REDDENING c C get index into color arrays IBVMS = (BVMSI +1 )/DBV + 1 IBVG = (BVGI +1 )/DBV + 1 C check if index is in range IF (IBVMS .LT. 1) IBVMS = 1 IF (IBVMS .GT. IBVMAX) IBVMS = IBVMAX IF (IBVG .LT. 1) IBVG = 1 IF (IBVG .GT. IBVMAX) IBVG = IBVMAX C accumulate component and total color arrays XBV(IBVMS,ID) = XBV(IBVMS,ID) + FMS XBV(IBVMS, 3) = XBV(IBVMS, 3) + FMS XBV(IBVG ,ID) = XBV(IBVG ,ID) + FG XBV(IBVG , 3) = XBV(IBVG , 3) + FG C C compute color means BVB( 3) = BVB( 3) + FMS*BVMSI + FG*BVGI BVB(ID) = BVB(ID) + FMS*BVMSI + FG*BVGI TBV( 3) = TBV( 3) + FG + FMS TBV(ID) = TBV(ID) + FG + FMS C compute fraction on giant branch XGB(ID) = XGB(ID) + FG XGB( 3) = XGB( 3) + FG RETURN END SUBROUTINE PRINT (IDS,XNM,DM,MLIM,XNMS,L2,B2,ABSM) IMPLICIT REAL (A-H,L,M,O-Z) REAL XNM(100), XNMS(100) C C print heading WRITE (8,1) 1 FORMAT ('0 m1',T11,'m2',T19,'m',T28,'A(m)', 1 T38,'N(