PROGRAM ABSDIM ! John Taylor 18/07/2004 ! Astrophysics Group, Keele University ! This inputs the observed quantities of an eclipsing binary and calcul- ! ates the absolute dimensions of the system and other useful quantities ! For each calculated quantity, the error arising from errors in each ! observed quantity is found and outputted, along with a total error. ! Total error is the individual errors added in quadrature. ! ! Written and checked: 20/7/2004 ! Improved array handling: 28/7/2004 ! Further distance indicators implemented: 4/8/2004 ! Tidied up: 22/11/04 ! Allowed facility of restricted use: 16/3/2005 ! Added Flower (1996) bolometric corrections: 20/8/2005 ! ! THE DATAFILES ARE ASSUMED TO BE IN DIRECTORY /home/jkt/tabs ! If they are in a different directory this must be changed in the code ! (just search for /home/jkt/ and replace the filenames). !======================================================================= implicit none integer NPAR ! Number of parameters input integer NRES ! Number of results calculated integer MH ! Metal abundance [M/H] flag integer WHAT ! Flag to say what to do parameter ( NPAR=29 , NRES=118 ) real*8 PAR(NPAR) ! Input parameters real*8 PARERR(NPAR) ! Uncertainties in parameters real*8 RES(NRES) ! Results calculated real*8 EHIGH(NRES),ELOW(NRES) ! Perturbed result values real*8 ERRORS(NPAR,NRES) ! Individual uncertainties real*8 SUM ! Summing variable real*8 RESERR(NRES) ! Overall result uncertainties character NAME(NRES)*10 ! Names of results calculated character UNIT(NRES)*5 ! Units the results are in real*8 FIT(NRES) ! Cosmic error fit uncertainty character BCTABLES(9)*35 ! Bolometric corr. table names integer NCODE,NBESSELL,NGIRARDI ! Numbers of BCs in tables real*8 CODE (100,3) ! BC data tabulated by Code real*8 BESSELL(1000,4) ! BC data tabulated by Bessell real*8 GIRARDI(1000,12) ! BC data tabulated by Girardi integer i,j,ERROR ! Loop counters and error flag common / BCDATA / CODE,BESSELL,GIRARDI,NCODE,NBESSELL,NGIRARDI BCTABLES(1) = "/home/jkt/tabs/absdim-CodeBCs.dat " BCTABLES(2) = "/home/jkt/tabs/absdim-BesselBCs.dat" BCTABLES(3) = "/home/jkt/tabs/absdim-GirBCs-25.dat" BCTABLES(4) = "/home/jkt/tabs/absdim-GirBCs-20.dat" BCTABLES(5) = "/home/jkt/tabs/absdim-GirBCs-15.dat" BCTABLES(6) = "/home/jkt/tabs/absdim-GirBCs-10.dat" BCTABLES(7) = "/home/jkt/tabs/absdim-GirBCs-05.dat" BCTABLES(8) = "/home/jkt/tabs/absdim-GirBCs+00.dat" BCTABLES(9) = "/home/jkt/tabs/absdim-GirBCs+05.dat" NAME( 1) = "Mass ratio" ; UNIT( 1) = " " ; FIT( 1) = 0.0d0 NAME( 2) = "M1sin(i)^3" ; UNIT( 2) = "Msun" ; FIT( 2) = 0.0d0 NAME( 3) = "M2sin(i)^3" ; UNIT( 3) = "Msun" ; FIT( 3) = 0.0d0 NAME( 4) = "a sin(i) " ; UNIT( 4) = "Rsun" ; FIT( 4) = 0.0d0 NAME( 5) = "Mass 1 " ; UNIT( 5) = "Msun" ; FIT( 5) = 0.0d0 NAME( 6) = "Mass 2 " ; UNIT( 6) = "Msun" ; FIT( 6) = 0.0d0 NAME( 7) = "Separation" ; UNIT( 7) = "Rsun" ; FIT( 7) = 0.0d0 NAME( 8) = "Radius 1 " ; UNIT( 8) = "Rsun" ; FIT( 8) = 0.0d0 NAME( 9) = "Radius 2 " ; UNIT( 9) = "Rsun" ; FIT( 9) = 0.0d0 NAME(10) = "log(g1) " ; UNIT(10) = "cm/s" ; FIT(10) = 0.0d0 NAME(11) = "log(g2) " ; UNIT(11) = "cm/s" ; FIT(11) = 0.0d0 NAME(12) = "V synch 1 " ; UNIT(12) = "km/s" ; FIT(12) = 0.0d0 NAME(13) = "V synch 2 " ; UNIT(13) = "km/s" ; FIT(13) = 0.0d0 NAME(14) = "time(sync)" ; UNIT(14) = "year" ; FIT(14) = 0.0d0 NAME(15) = "time(circ)" ; UNIT(15) = "year" ; FIT(15) = 0.0d0 NAME(16) = "log(L) 1 " ; UNIT(16) = "Lsun" ; FIT(16) = 0.0d0 NAME(17) = "log(L) 2 " ; UNIT(17) = "Lsun" ; FIT(17) = 0.0d0 NAME(18) = "M(bol) 1 " ; UNIT(18) = "mag " ; FIT(18) = 0.0d0 NAME(19) = "M(bol) 2 " ; UNIT(19) = "mag " ; FIT(19) = 0.0d0 NAME(20) = "M(V) 1 " ; UNIT(20) = "mag " ; FIT(20) = 0.0d0 NAME(21) = "M(V) 2 " ; UNIT(21) = "mag " ; FIT(21) = 0.0d0 NAME(22) = "M(V) total" ; UNIT(22) = "mag " ; FIT(22) = 0.0d0 NAME(23) = "distance 1" ; UNIT(23) = "pc " ; FIT(23) = 0.0d0 NAME(24) = "distance 2" ; UNIT(24) = "pc " ; FIT(24) = 0.0d0 NAME(25) = "d(overall)" ; UNIT(25) = "pc " ; FIT(25) = 0.0d0 NAME(26) = "M(K) 1 " ; UNIT(26) = "mag " ; FIT(26) = 0.0d0 NAME(27) = "M(K) 2 " ; UNIT(27) = "mag " ; FIT(27) = 0.0d0 NAME(28) = "M(K) total" ; UNIT(28) = "mag " ; FIT(28) = 0.0d0 NAME(29) = "distance 1" ; UNIT(29) = "pc " ; FIT(29) = 0.0d0 NAME(30) = "distance 2" ; UNIT(30) = "pc " ; FIT(30) = 0.0d0 NAME(31) = "d(overall)" ; UNIT(31) = "pc " ; FIT(31) = 0.0d0 NAME(32) = "M(U) 1 " ; UNIT(32) = "mag " ; FIT(32) = 0.0d0 NAME(33) = "M(U) 2 " ; UNIT(33) = "mag " ; FIT(33) = 0.0d0 NAME(34) = "M(U) total" ; UNIT(34) = "mag " ; FIT(34) = 0.0d0 NAME(35) = "distance 1" ; UNIT(35) = "pc " ; FIT(35) = 0.0d0 NAME(36) = "distance 2" ; UNIT(36) = "pc " ; FIT(36) = 0.0d0 NAME(37) = "d(overall)" ; UNIT(37) = "pc " ; FIT(37) = 0.0d0 NAME(38) = "M(B) 1 " ; UNIT(38) = "mag " ; FIT(38) = 0.0d0 NAME(39) = "M(B) 2 " ; UNIT(39) = "mag " ; FIT(39) = 0.0d0 NAME(40) = "M(B) total" ; UNIT(40) = "mag " ; FIT(40) = 0.0d0 NAME(41) = "distance 1" ; UNIT(41) = "pc " ; FIT(41) = 0.0d0 NAME(42) = "distance 2" ; UNIT(42) = "pc " ; FIT(42) = 0.0d0 NAME(43) = "d(overall)" ; UNIT(43) = "pc " ; FIT(43) = 0.0d0 NAME(44) = "M(V) 1 " ; UNIT(44) = "mag " ; FIT(44) = 0.0d0 NAME(45) = "M(V) 2 " ; UNIT(45) = "mag " ; FIT(45) = 0.0d0 NAME(46) = "M(V) total" ; UNIT(46) = "mag " ; FIT(46) = 0.0d0 NAME(47) = "distance 1" ; UNIT(47) = "pc " ; FIT(47) = 0.0d0 NAME(48) = "distance 2" ; UNIT(48) = "pc " ; FIT(48) = 0.0d0 NAME(49) = "d(overall)" ; UNIT(49) = "pc " ; FIT(49) = 0.0d0 NAME(50) = "M(R) 1 " ; UNIT(50) = "mag " ; FIT(50) = 0.0d0 NAME(51) = "M(R) 2 " ; UNIT(51) = "mag " ; FIT(51) = 0.0d0 NAME(52) = "M(R) total" ; UNIT(52) = "mag " ; FIT(52) = 0.0d0 NAME(53) = "distance 1" ; UNIT(53) = "pc " ; FIT(53) = 0.0d0 NAME(54) = "distance 2" ; UNIT(54) = "pc " ; FIT(54) = 0.0d0 NAME(55) = "d(overall)" ; UNIT(55) = "pc " ; FIT(55) = 0.0d0 NAME(56) = "M(I) 1 " ; UNIT(56) = "mag " ; FIT(56) = 0.0d0 NAME(57) = "M(I) 2 " ; UNIT(57) = "mag " ; FIT(57) = 0.0d0 NAME(58) = "M(I) total" ; UNIT(58) = "mag " ; FIT(58) = 0.0d0 NAME(59) = "distance 1" ; UNIT(59) = "pc " ; FIT(59) = 0.0d0 NAME(60) = "distance 2" ; UNIT(60) = "pc " ; FIT(60) = 0.0d0 NAME(61) = "d(overall)" ; UNIT(61) = "pc " ; FIT(61) = 0.0d0 NAME(62) = "M(J) 1 " ; UNIT(62) = "mag " ; FIT(62) = 0.0d0 NAME(63) = "M(J) 2 " ; UNIT(63) = "mag " ; FIT(63) = 0.0d0 NAME(64) = "M(J) total" ; UNIT(64) = "mag " ; FIT(64) = 0.0d0 NAME(65) = "distance 1" ; UNIT(65) = "pc " ; FIT(65) = 0.0d0 NAME(66) = "distance 2" ; UNIT(66) = "pc " ; FIT(66) = 0.0d0 NAME(67) = "d(overall)" ; UNIT(67) = "pc " ; FIT(67) = 0.0d0 NAME(68) = "M(H) 1 " ; UNIT(68) = "mag " ; FIT(68) = 0.0d0 NAME(69) = "M(H) 2 " ; UNIT(69) = "mag " ; FIT(69) = 0.0d0 NAME(70) = "M(H) total" ; UNIT(70) = "mag " ; FIT(70) = 0.0d0 NAME(71) = "distance 1" ; UNIT(71) = "pc " ; FIT(71) = 0.0d0 NAME(72) = "distance 2" ; UNIT(72) = "pc " ; FIT(72) = 0.0d0 NAME(73) = "d(overall)" ; UNIT(73) = "pc " ; FIT(73) = 0.0d0 NAME(74) = "M(K) 1 " ; UNIT(74) = "mag " ; FIT(74) = 0.0d0 NAME(75) = "M(K) 2 " ; UNIT(75) = "mag " ; FIT(75) = 0.0d0 NAME(76) = "M(K) total" ; UNIT(76) = "mag " ; FIT(76) = 0.0d0 NAME(77) = "distance 1" ; UNIT(77) = "pc " ; FIT(77) = 0.0d0 NAME(78) = "distance 2" ; UNIT(78) = "pc " ; FIT(78) = 0.0d0 NAME(79) = "d(overall)" ; UNIT(79) = "pc " ; FIT(79) = 0.0d0 NAME(80) = "M(V) 1 " ; UNIT(80) = "mag " ; FIT(80) = 0.0d0 NAME(81) = "M(V) 2 " ; UNIT(81) = "mag " ; FIT(81) = 0.0d0 NAME(82) = "M(V) total" ; UNIT(82) = "mag " ; FIT(82) = 0.0d0 NAME(83) = "distance 1" ; UNIT(83) = "pc " ; FIT(83) =0.025d0 NAME(84) = "distance 2" ; UNIT(84) = "pc " ; FIT(84) =0.025d0 NAME(85) = "d(overall)" ; UNIT(85) = "pc " ; FIT(85) =0.025d0 NAME(86) = "M(V) 1 " ; UNIT(86) = "mag " ; FIT(86) = 0.0d0 NAME(87) = "M(V) 2 " ; UNIT(87) = "mag " ; FIT(87) = 0.0d0 NAME(88) = "M(V) total" ; UNIT(88) = "mag " ; FIT(88) = 0.0d0 NAME(89) = "distance 1" ; UNIT(89) = "pc " ; FIT(89) =0.025d0 NAME(90) = "distance 2" ; UNIT(90) = "pc " ; FIT(90) =0.025d0 NAME(91) = "d(overall)" ; UNIT(91) = "pc " ; FIT(91) =0.025d0 NAME(92) = "d(U) 1 " ; UNIT(92) = "pc " ; FIT(92) =0.059d0 NAME(93) = "d(U) 2 " ; UNIT(93) = "pc " ; FIT(93) =0.059d0 NAME(94) = "d(U) total" ; UNIT(94) = "pc " ; FIT(94) =0.059d0 NAME(95) = "d(B) 1 " ; UNIT(95) = "pc " ; FIT(95) =0.063d0 NAME(96) = "d(B) 2 " ; UNIT(96) = "pc " ; FIT(96) =0.063d0 NAME(97) = "d(B) total" ; UNIT(97) = "pc " ; FIT(97) =0.063d0 NAME(98) = "d(V) 1 " ; UNIT(98) = "pc " ; FIT(98) =0.059d0 NAME(99) = "d(V) 2 " ; UNIT(99) = "pc " ; FIT(99) =0.059d0 NAME(100)= "d(V) total" ; UNIT(100)= "pc " ; FIT(100)=0.059d0 NAME(101)= "d(R) 1 " ; UNIT(101)= "pc " ; FIT(101)=0.048d0 NAME(102)= "d(R) 2 " ; UNIT(102)= "pc " ; FIT(102)=0.048d0 NAME(103)= "d(R) total" ; UNIT(103)= "pc " ; FIT(103)=0.048d0 NAME(104)= "d(I) 1 " ; UNIT(104)= "pc " ; FIT(104)=0.023d0 NAME(105)= "d(I) 2 " ; UNIT(105)= "pc " ; FIT(105)=0.023d0 NAME(106)= "d(I) total" ; UNIT(106)= "pc " ; FIT(106)=0.023d0 NAME(107)= "d(J) 1 " ; UNIT(107)= "pc " ; FIT(107)=0.012d0 NAME(108)= "d(J) 2 " ; UNIT(108)= "pc " ; FIT(108)=0.012d0 NAME(109)= "d(J) total" ; UNIT(109)= "pc " ; FIT(109)=0.012d0 NAME(110)= "d(H) 1 " ; UNIT(110)= "pc " ; FIT(110)=0.014d0 NAME(111)= "d(H) 2 " ; UNIT(111)= "pc " ; FIT(111)=0.014d0 NAME(112)= "d(H) total" ; UNIT(112)= "pc " ; FIT(112)=0.014d0 NAME(113)= "d(K) 1 " ; UNIT(113)= "pc " ; FIT(113)= 0.01d0 NAME(114)= "d(K) 2 " ; UNIT(114)= "pc " ; FIT(114)= 0.01d0 NAME(115)= "d(K) total" ; UNIT(115)= "pc " ; FIT(115)= 0.01d0 NAME(116)= "d(L) 1 " ; UNIT(116)= "pc " ; FIT(116)= 0.01d0 NAME(117)= "d(L) 2 " ; UNIT(117)= "pc " ; FIT(117)= 0.01d0 NAME(118)= "d(L) total" ; UNIT(118)= "pc " ; FIT(118)= 0.01d0 do i = 1,NPAR PAR(i) = -100.0d0 end do do i = 1,NRES RES(i) = -100.0d0 end do do i = 1,NPAR do j = 1,NRES ERRORS(i,j) = 0.0d0 end do end do WHAT = 0 !----------------------------------------------------------------------- write (*,*) " " write (*,'(A40,A40)') "ABSDIM (John Taylor 16/7/04) This reads ", & "in spectroscopic and photometric results" write (*,'(A40,A40)') "for an eclipsing binary and outputs its ", & "absolute dimensions, distance and other." write (*,'(A40,A40)') "Enter the orbital period, K1, K2, orbita", & "l eccentricity, inclination, r1 and r2: " read (*,*) (PAR(i),i=1,7) write (*,'(A40,A40)') "Enter the absolute uncertainties of the ", & "above quantities, in the same order: " read (*,*) (PARERR(i),i=1,7) write (*,'(A40,A40)') "Enter the two Teffs, reddening E(B-V) an", & "d [M/H] (from 1 for -2.5 to 7 for +0.5)" write (*,'(A40,A40)') "(if no photometric results are required ", & "then put Teff1 or Teff2 less than zero):" read (*,*) (PAR(i),i=8,10),MH if ( PAR(8) < 0.0d0 .or. PAR(9) < 0.0d0 ) GOTO 200 write (*,'(A40,A40)') "Enter the absolute uncertainties in the ", & "quantities Teff1, Teff2 and reddening: " read (*,*) (PARERR(i),i=8,10) write (*,'(A40,A40)') "Enter apparent magnitudes of system in U", & "BVRIJHKL (-100.0 for unavailable data): " read (*,*) (PAR(i),i=11,19) write (*,'(A40,A40)') "Enter the absolute uncertainties of the ", & "above quantities, in the same order: " read (*,*) (PARERR(i),i=11,19) write (*,'(A40,A40)') "Enter the light ratios of the system in ", & "UBVRIJHKL (-100.0 for unavailable data):" read (*,*) (PAR(i),i=20,28) write (*,'(A40,A40)') "Enter the absolute uncertainties of the ", & "above quantities, in the same order: " read (*,*) (PARERR(i),i=20,28) WHAT = 1 !----------------------------------------------------------------------- open (56,file=BCTABLES(1),status="old",iostat=ERROR) if (ERROR /= 0) write (*,'(A35,A35)') & "### ERROR: unable to open the file ",BCTABLES(1) do i = 1,100 read (56,*,iostat=ERROR) (CODE(i,j),j=1,3) if ( ERROR /= 0 ) exit end do NCODE = i - 1 close (unit=56) open (56,file=BCTABLES(2),status="old",iostat=ERROR) if (ERROR /= 0) write (*,'(A35,A35)') & "### ERROR: unable to open the file ",BCTABLES(2) do i = 1,1000 read (56,*,iostat=ERROR) (BESSELL(i,j),j=1,4) if ( ERROR /= 0 ) exit end do NBESSELL = i - 1 close (unit=56) open (unit=56,file=BCTABLES(MH+2),status="old",iostat=ERROR) if ( ERROR /= 0 ) write(*,'(A35,A35)') & "### ERROR: unable to open the file ",BCTABLES(MH+2) do i = 1,1000 read (56,*,iostat=ERROR) (GIRARDI(i,j),j=1,12) if ( ERROR /= 0 ) exit end do NGIRARDI = i - 1 close (unit=56) !----------------------------------------------------------------------- 200 write (*,*) " " write(*,'(A43,A64)')"Quantity result error ", &" dP dK1 dK2 de di dr1 dr2 dT1 dT2 EBV mag dLR fit" write(*,'(A43,A64)')"-------------------------------------------", &"----------------------------------------------------------------" ! First calculate the results for the given parameters. Then ! for each parameter, perturb its value by plus and by minus ! its stated uncertainty to find the uncertainty in each ! result due only to the uncertainty in that parameter. CALL CALC (PAR,NPAR,NRES,RES) do i = 1,NPAR if ( WHAT == 0 .and. i >= 16 ) exit if ( PAR(i) > -99.0d0 ) then PAR(i) = PAR(i) + PARERR(i) CALL CALC (PAR,NPAR,NRES,EHIGH) PAR(i) = PAR(i) - (2.0d0 * PARERR(i)) CALL CALC (PAR,NPAR,NRES,ELOW) PAR(i) = PAR(i) + PARERR(i) do j = 1,NRES ERRORS(i,j) = 0.5d0 * abs(EHIGH(j) - ELOW(j)) end do end if end do ! Now the individual errors need to be sorted for output. ! First, the fractional error due to cosmic scatter in any ! calibrations is deposited into ERRORS(28). ! Second, the total error in each result is found by adding ! in quadrature the individual errors found by individually ! perturbing every input parameter by its error. ! Third, the individual errors for a result due to each ! parameter are then scaled to be a fraction of the total ! error found for that result. do i = 1,NRES ERRORS(29,i) = FIT(i) * RES(i) end do do j = 1,NRES SUM = 0.0d0 do i = 1,NPAR SUM = SUM + ERRORS(i,j)**2 end do RESERR(j) = sqrt(SUM) do i = 1,NPAR ERRORS(i,j) = ERRORS(i,j) / RESERR(j) end do end do ! As only one magnitude or flux ratio has bearing on each ! result, the uncertainties for only the relevant one need ! to be outputted each time. Achieve this by putting the ! sum of the individual errors from parameters 11 to 19 ! (UBVRIJHKL magnitudes) into ERRORS(11,i), and likewise ! for 20 to 27 going into 12. Then 28 goes into 13. do i = 1,NRES ERRORS(11,i) = sqrt(ERRORS(11,i)**2 + ERRORS(12,i)**2 + & ERRORS(13,i)**2 + ERRORS(14,i)**2 + & ERRORS(15,i)**2 + ERRORS(16,i)**2 + & ERRORS(17,i)**2 + ERRORS(18,i)**2 + ERRORS(19,i)**2) ERRORS(12,i) = sqrt(ERRORS(20,i)**2 + ERRORS(21,i)**2 + & ERRORS(22,i)**2 + ERRORS(23,i)**2 + & ERRORS(24,i)**2 + ERRORS(25,i)**2 + & ERRORS(26,i)**2 + ERRORS(27,i)**2 + ERRORS(28,i)**2) ERRORS(13,i) = ERRORS(29,i) end do ! Now output the results, their overall uncertainties, units ! and the individual contributions from each parameter. do j = 1,NRES if ( WHAT == 0 .and. j >= 16 ) exit if ( RES(j) > -99.0d0 ) then write (*,500) NAME(j),RES(j),RESERR(j),UNIT(j) do i = 1,13 if ( ERRORS(i,j) >= 0.9995d0 ) then write (*,501) ERRORS(i,j) else if ( ERRORS(i,j) > 0.0005d0 ) then write (*,502) ERRORS(i,j) else write (*,503) " ---" end if end do write (*,*) " " ! Line break end if if ( j == 19 ) write (*,'(A52,A52)') & "Absolute magnitudes and distances derived using bolo", & "metric corrections of Bessell (1998A+A...333..231B):" if ( j == 31 ) write (*,'(A52,A52)') & "Absolute magnitudes and distances derived using bolo", & "metric corrections of Girardi (2002A+A...391..195G):" if ( j == 79 ) write (*,'(A52,A52)') & "Absolute magnitudes and distances using empirical bo", & "lometric corrections of Code (1976ApJ...203..417C): " if ( j == 86 ) write (*,'(A52,A52)') & "Absolute mag and distances with empirical bolometric", & " correction formula of Flower (1996ApJ...469..355F):" if ( j == 92 ) write (*,'(A52,A52)') & "Distances using the UBVRIJHKL surface brightness-Tef", & "f relations of Kervella et al (2004A+A...426..297K):" if ( j == 92 ) then if (PAR(8)>1.0d4.or.PAR(8)<3.6d3) write(*,'(A52,A35,F8.1,A7)') & "### Warning: the Kervella relations are valid for 36", & "00 < Teff < 10000 K, but Teff(A) is",PAR(8)," Kelvin" if (PAR(9)>1.0d4.or.PAR(9)<3.6d3) write(*,'(A52,A35,F8.1,A7)') & "### Warning: the Kervella relations are valid for 36", & "00 < Teff < 10000 K, but Teff(B) is",PAR(9)," Kelvin" end if end do write (*,*) " " 500 FORMAT (A10," =",F11.6," +/-",F10.6,1X,A4,$) 501 FORMAT (1X,F4.2,$) 502 FORMAT (1X,F4.3,$) 503 FORMAT (A5,$) END PROGRAM ABSDIM !======================================================================= SUBROUTINE CALC (IN,NIN,NOUT,OUT) ! IN: ! (1) orbital period (days) ! (2),(3) velocity semiamplitudes of the stars (kms/s) ! (4),(5) orbital eccentricity and inclination (degrees) ! (6),(7) fractional radii of the stars ! (8),(9) effective temperatures of the stars (Kelvin) ! (10) reddening of the system E(B-V) (magnitudes) ! (11)-(19) apparent system magnitudes in U B V R I J H K L ! (20)-(28) flux ratios in U B V R I J H K L ! (29) cosmic scatter in calibrations ! OUT: ! (1) the mass ratio ! (2),(3) minimum masses of the stars (Msun) ! (4) the projected separation of the stars (Rsun) ! (5),(6) absolute mass of the stars (Msun) ! (7) absolute separation of the star centres (Rsun) ! (8),(9) absolute radii of the stars (Rsun) ! (10),(11) surface gravities, log(g), in log(cm/s) ! (12),(13) synchronous rotational velocities (km/s) ! (14) synchronisation timescale (Zahn approx), log(yr) ! (15) circularisation timescale (Zahn approx), log(yr) ! (16),(17) luminosities of the stars (W) ! (18),(19) absolute bolometric magnitudes of the stars ! (20)-(25) Bessell BCs: MV1,MV2,MVtot,DIST1,DIST2,DISTtot ! (26)-(31) Bessell BCs: MK1,MK2,MKtot,DIST1,DIST2,DISTtot ! (32)-(79) 8 lots of (20)-(25), for Girardi BCs in UBVRIJHK ! (80)-(85) One lot of (20)-(25) for Code empirical BCs in V ! (86)-(91) One lot of (20)-(25) for Flower (1996) V BCs ! (92)-(118) Distances from Kervella SB-Teff calibrations implicit none integer NIN,NOUT ! IN: size of in and out array real*8 IN(NIN) ! IN: the inputs parameters real*8 OUT(NOUT) ! OUT: the calculated results real*8 PI,DEG2RAD,SB ! LOCAL: constants real*8 AU,CONST1,CONST2 ! LOCAL: constants real*8 SINI,EFAC,KSUM ! LOCAL: calculation quantitis real*8 MSUN,RSUN,LSUN,MBOLSUN ! LOCAL: properties of the Sun real*8 BC1,BC2 ! LOCAL: bolometric correction real*8 MAG(9),MAG1(9),MAG2(9) ! LOCAL: magnitudes (A+B,A,B) integer BCCOL(8) ! LOCAL: columns to pick integer i,j ! LOCAL: loop counters data BCCOL/4,5,7,8,9,10,11,12/ PI = atan(1.0d0) * 4.0d0 ! pi to the machine precision DEG2RAD = PI / 180.0d0 ! degrees-radians conversion SB = 5.670400d-8 CONST1 = 1.036149d-7 ! constant in mass formula CONST2 = 1.9758d-2 ! constant in a*sin(i) formula AU = 1.495979d11 ! Astronomical Unit in metres MSUN = 1.9891d26 ! The mass (kg), radius (m), RSUN = 6.9599d8 ! luminosity (W) and absolute LSUN = 3.826d26 ! bolometric magnitude of the MBOLSUN = 4.75 ! Sun have come from Zombeck. SINI = sin( IN(5) * DEG2RAD ) ! sine(inclination) EFAC = 1.0d0 - IN(4)**2 ! eccentricity factor KSUM = IN(2) + IN(3) ! velocity semiamplitude sum !----------------------------------------------------------------------- OUT(1) = IN(2) / IN(3) ! Mass ratio OUT(2) = CONST1 * EFAC**1.5d0 * KSUM**2 * IN(3) * IN(1) OUT(3) = CONST1 * EFAC**1.5d0 * KSUM**2 * IN(2) * IN(1) OUT(4) = CONST2 * sqrt(EFAC) * KSUM * IN(1) ! Projected sep'n OUT(5) = OUT(2) / SINI**3 ! Mass of star 1 OUT(6) = OUT(3) / SINI**3 ! Mass of star 2 OUT(7) = OUT(4) / SINI ! Separation if ( IN(7) < 0.0d0 ) then OUT(8) = OUT(7) * IN(6) / (1.0d0 + abs(IN(7))) OUT(9) = OUT(7) * IN(6) * abs(IN(7)) / (1.0d0 + abs(IN(7))) else OUT(8) = IN(6) * OUT(7) ! Radius of star 1 OUT(9) = IN(7) * OUT(7) ! Radius of star 2 end if OUT(10) = 4.438d0 + log10(OUT(5)) - (2.0d0 * log10(OUT(8))) OUT(11) = 4.438d0 + log10(OUT(6)) - (2.0d0 * log10(OUT(9))) OUT(12) = 50.58d0 * OUT(8) / IN(1) OUT(13) = 50.58d0 * OUT(9) / IN(1) ! Calculate the timescales of rotational synchronisation and ! orbital circularisation using Zahn's approximations for ! stars with convective envelopes (1977A+A....57..383). OUT(14) = 1.0d4 * ((1.0d0+OUT(1))/(2.0d0*OUT(1)))**2 * IN(1)**4 OUT(14) = log10( OUT(14) ) OUT(15) = 1.0d6 / OUT(1) * ((1.0d0*OUT(1))/2.0d0)**(5.0d0/3.0d0) OUT(15) = log10( OUT(15) * IN(1)**(16.0d0/3.0d0) ) ! Calculate the luminosities and absolute bolometric magnit- ! udes using traditional data for the Sun (Zombeck 1990). OUT(16) = log10((4.0d0*PI*SB * (OUT(8)*RSUN)**2 * IN(8)**4)/LSUN) OUT(17) = log10((4.0d0*PI*SB * (OUT(9)*RSUN)**2 * IN(9)**4)/LSUN) OUT(18) = (-2.5d0 * OUT(16)) + MBOLSUN OUT(19) = (-2.5d0 * OUT(17)) + MBOLSUN !----------------------------------------------------------------------- CALL DEREDDEN (IN(10),IN(11),IN(20),MAG,MAG1,MAG2) ! Calculate the distance using the bolometric corrections of ! Bessell, Castelli & Plez (1998A+A...333..231B) for V and K ! Must adopt LSUN=3.855e26 and MBOLSUN=4.74 for consistency. ! This requires recalculation of L and MBOL for the stars. if ( MAG(3) > -99.0d0 ) then CALL BCDISTANCE (1,IN(8),IN(9),OUT(10),OUT(11),OUT(8),OUT(9), & 3.855d26,4.74d0,MAG(3),MAG1(3),MAG2(3),1,2,4, & OUT(20),OUT(21),OUT(22),OUT(23),OUT(24),OUT(25)) end if if ( MAG(8) > -99.0d0 ) then CALL BCDISTANCE (1,IN(8),IN(9),OUT(10),OUT(11),OUT(8),OUT(9), & 3.855d26,4.74d0,MAG(8),MAG1(8),MAG2(8),1,2,3, & OUT(26),OUT(27),OUT(28),OUT(29),OUT(30),OUT(31)) end if ! Calculate the distance using the bolometric corrections of ! Girardi et al (2002A+A...391..195G) for Ux,Bx,V,R,I,J,H,K. ! Must adopt LSUN=3.844e26 and MBOLSUN=4.77 for consistency. do i = 1,8 j = 26 + 6*i if ( MAG(i) > -99.0d0 ) then CALL BCDISTANCE (2,IN(8),IN(9),OUT(10),OUT(11),OUT(8),OUT(9), & 3.844d26,4.77d0,MAG(i),MAG1(i),MAG2(i),2,3,BCCOL(i), & OUT(j),OUT(j+1),OUT(j+2),OUT(j+3),OUT(j+4),OUT(j+5)) end if end do ! Calculate the distance using the empirical bolometric cor- ! rections of Code et al (1976ApJ...203..417C) for V filter. ! The solar data is the same as adopted above. if ( MAG(3) > -99.0d0 ) then CALL BCDISTANCE (3,IN(8),IN(9),OUT(10),OUT(11),OUT(8),OUT(9), & LSUN,MBOLSUN,MAG(3),MAG1(3),MAG2(3),1,2,3, & OUT(80),OUT(81),OUT(82),OUT(83),OUT(84),OUT(85)) end if ! Calculate the distance using the empirical bolometric co- ! rection formula of Flower (1996ApJ...469..355F) for V band ! Adopt the solar data as above, as Flower doesn't say. if ( MAG(3) > -99.0d0 ) then CALL FLOWER (IN(8),IN(9),OUT(8),OUT(9),LSUN,MBOLSUN,MAG(3),MAG1( & 3),MAG2(3),OUT(86),OUT(87),OUT(88),OUT(89),OUT(90),OUT(91)) end if !----------------------------------------------------------------------- ! Calculate the distance using surface brightness-effective ! temperature relations of Kervella (astro-ph/0404180) do i = 1,9 CALL KERV (IN(8),IN(9),OUT(8),OUT(9),i,MAG(i),MAG1(i),MAG2(i), & RSUN,AU,OUT(89+3*i),OUT(90+3*i),OUT(91+3*i)) end do !----------------------------------------------------------------------- END SUBROUTINE CALC !======================================================================= SUBROUTINE KERV (TEFF1,TEFF2,RADIUS1,RADIUS2,WHICH,MAG,MAG1,MAG2, & RSUN,AU,DIST1,DIST2,DIST) ! Calculates the distance using surface brightness-effective ! temperature relations of Kervella (astro-ph/0404180) for ! the broad-band UBVRIJHKL filters. The distance to each star ! and to the whole system is calculated. implicit none real*8 TEFF1,TEFF2,RADIUS1,RADIUS2 ! IN: Teffs and radii of stars integer WHICH ! IN: 1 to 9 for UBVRIJHKL real*8 MAG,MAG1,MAG2 ! IN: app mag of dEB, stars real*8 RSUN,AU ! IN: solar constants real*8 DIST1,DIST2,DIST ! OUT: distances to stars, dEB real*8 LOGTHETA1,LOGTHETA2 ! LOCAL: angular diameters real*8 ZMLD1,ZMLD2 ! LOCAL: surface brightnesses real*8 GETLOGZMLD ! LOCAL: FUNCTION real*8 FACTOR,HELPER ! LOCAL: helper variables if ( MAG1>=-99.0d0 .and. TEFF1 < 1.5d4 ) then LOGTHETA1 = GETLOGZMLD(TEFF1,WHICH) - 0.2d0*MAG1 DIST1 = 9.3048d0 * RADIUS1 / 10.0d0**LOGTHETA1 else DIST1 = -100.0d0 end if if ( MAG2>=-99.0d0 .and. TEFF2 < 1.5d4 ) then LOGTHETA2 = GETLOGZMLD(TEFF2,WHICH) - 0.2d0*MAG2 DIST2 = 9.3048d0 * RADIUS2 / 10.0d0**LOGTHETA2 else DIST2 = -100.0d0 end if if ( MAG >=-99.0d0 .and. TEFF1 < 1.5d4 .and. TEFF2 < 1.5d4 ) then FACTOR = 1000.0d0 * RSUN / AU ZMLD1 = 10.0d0 ** GETLOGZMLD(TEFF1,WHICH) ZMLD2 = 10.0d0 ** GETLOGZMLD(TEFF2,WHICH) HELPER = (RADIUS1*FACTOR/ZMLD1)**2 + (RADIUS2*FACTOR/ZMLD2)**2 DIST = 10.0d0**(0.2d0*MAG) * 2.0d0 * sqrt(HELPER) else DIST = -100.0d0 end if END SUBROUTINE KERV !======================================================================= DOUBLEPRECISION FUNCTION GETLOGZMLD (TEFF,WHICH) ! Given an effective temperature, calculates zero-magnitude ! angular diameter (mas) using the calibration of Kervella. implicit none real*8 TEFF,A(9),B(9),C(9) ! IN: Teff LOCAL: fit coefficients integer WHICH ! IN: Which filter (UBVRIJHKL) data A /5.6391d0,3.6753d0,3.0415d0,2.1394d0,0.9847d0, & 0.9598d0,1.1684d0,0.8470d0,0.6662d0/ data B /-46.4505d0,-30.9671d0,-25.4696d0,-18.0221d0,-8.7985d0, & -8.3451d0,-9.6156d0,-7.0790d0,-5.6609/ data C /96.0513d0,65.5421d0,53.7010d0,38.3497d0,19.9281d0, & 18.5204d0,20.2779d0,15.2731d0,12.4902/ GETLOGZMLD = A(WHICH)*log10(TEFF)**2+B(WHICH)*log10(TEFF)+C(WHICH) END FUNCTION GETLOGZMLD !======================================================================= !======================================================================= SUBROUTINE DEREDDEN (EBV,DEBMAG,FLUXRAT,MAG,MAG1,MAG2) ! Calculate the dereddened apparent magnitudes of the system ! and of each star (-10.0 means that no data exists). implicit none real*8 EBV ! IN: reddening E(B-V) real*8 DEBMAG(9) ! IN: system UBVRIJHKL apparent mags real*8 FLUXRAT(9) ! IN: flux ratios in UBVRIJHKL bands real*8 MAG(9),MAG1(9),MAG2(9) ! OUT: dereddened mags for dEB,1,2 real*8 AV(9) ! LOCAL: extinction ratios integer i ! LOCAL: loop counter data AV/5.0d0,4.2d0,3.2d0,2.5d0,2.0d0,1.0d0,0.6d0,0.5d0,0.4d0/ do i = 1,9 MAG(i) = -100.0d0 MAG1(i) = -100.0d0 MAG2(i) = -100.0d0 if ( DEBMAG(i) > -99.0d0 ) then MAG(i) = DEBMAG(i) - AV(i)*EBV if ( FLUXRAT(i) > -99.0d0 ) then MAG1(i) = MAG(i) + 2.5d0 * log10(1.0d0+FLUXRAT(i)) MAG2(i) = MAG(i)+2.5d0*log10((1.0d0+FLUXRAT(i))/FLUXRAT(i)) end if end if end do END SUBROUTINE DEREDDEN !======================================================================= SUBROUTINE BCDISTANCE (BCTAB,TEFF1,TEFF2,LOGG1,LOGG2,RADIUS1, & RADIUS2,LSUN,MBOLSUN,MAG,MAG1,MAG2,CTEFF, & CLOGG,CBC,MV1,MV2,MVTOT,DIST1,DIST2,DIST) implicit none real*8 TEFF1,TEFF2,LOGG1,LOGG2 ! IN: stellar parameters real*8 RADIUS1,RADIUS2 ! IN: stellar radii in Rsun real*8 LSUN,MBOLSUN ! IN: solar parametrs to adopt real*8 MAG,MAG1,MAG2 ! IN: mag of dEB, star1, star2 integer CTEFF,CLOGG,CBC ! IN: column numbers in file real*8 MV1,MV2,MVTOT ! OUT: M(V) for 1,2,dEB real*8 DIST1,DIST2,DIST ! OUT: distances to 1,2,dEB real*8 LOGL1,LOGL2,MBOL1,MBOL2 ! LOCAL: radiative parameters real*8 BC1,BC2 ! LOCAL: bolometric corections real*8 PI,SB,RSUN ! LOCAL: astro constants integer BCTAB PI = atan(1.0d0) * 4.0d0 ! Pi to the machine precision SB = 5.670400d-8 ! Stefan-Boltzmann constant RSUN = 6.9599d8 ! Radius of Sun from Zombeck DIST1 = -100.0d0 ! Initialise DIST2 = -100.0d0 ! Initialise DIST = -100.0d0 ! Initialise LOGL1 = log10( (4.0d0*PI*SB*(RADIUS1*RSUN)**2*TEFF1**4) / LSUN ) LOGL2 = log10( (4.0d0*PI*SB*(RADIUS2*RSUN)**2*TEFF2**4) / LSUN ) MBOL1 = (-2.5d0 * LOGL1) + MBOLSUN MBOL2 = (-2.5d0 * LOGL2) + MBOLSUN CALL BCGET (BCTAB,TEFF1,LOGG1,CTEFF,CLOGG,CBC,BC1) CALL BCGET (BCTAB,TEFF2,LOGG2,CTEFF,CLOGG,CBC,BC2) MV1 = MBOL1 - BC1 MV2 = MBOL2 - BC2 MVTOT = -2.5d0 * log10(10.0d0**(-0.4d0*MV1)+10.0d0**(-0.4d0*MV2)) if ( MAG1 > -99.0d0 ) DIST1 = 10.0d0 ** ( (MAG1-MV1+5.0d0)/5.0d0 ) if ( MAG2 > -99.0d0 ) DIST2 = 10.0d0 ** ( (MAG2-MV2+5.0d0)/5.0d0 ) if ( MAG > -99.0d0 ) DIST = 10.0d0 ** ( (MAG-MVTOT+5.0d0)/5.0d0 ) END SUBROUTINE BCDISTANCE !======================================================================= SUBROUTINE BCGET (BCTAB,TEFF,LOGG,CTEFF,CLOGG,CBC,BC) ! Given an effective temperature and a surface gravity, read ! ina file containing bolometric corrections and interpolate ! to find an appropriate correction for the Teff and log(g). ! Interpolation is linear and two-dimensional in Teff,log(g) implicit none integer BCTAB ! IN: 1 Bessell, 2 Girardi, 3 Code real*8 TEFF,LOGG ! IN: stellar parameters integer CTEFF,CLOGG,CBC ! IN: columns in file to use real*8 BC ! OUT: bolometric correction real*8 LOWTEFF,HIGHTEFF ! LOCAL: values which bracket Teff real*8 LOWLOGG,HIGHLOGG ! LOCAL: values which bracket log(g) real*8 INTDATA(4) ! LOCAL: data for interpolation real*8 WEIGHTS(4) ! LOCAL: weights for interpolation integer i,j ! LOCAL: loop counters real*8 CODE(100,3) ! BC data tabulated by Code et al real*8 BESSELL(1000,4) ! BC data tabulated by Bessell real*8 GIRARDI(1000,12) ! BC data tabulated by Girardi integer NCODE,NBESSELL,NGIRARDI,NUMDATA real*8 DATA(1000,12) common / BCDATA / CODE,BESSELL,GIRARDI,NCODE,NBESSELL,NGIRARDI data INTDATA / -10.0d0,-10.0d0,-10.0d0,-10.0d0 / LOWTEFF = 0.0d0 ; LOWLOGG = 0.0d0 HIGHTEFF = 1.0d8 ; HIGHLOGG = 10.0d0 if ( BCTAB == 1 ) then do i = 1,NBESSELL do j = 1,4 DATA(i,j) = BESSELL(i,j) end do end do NUMDATA = NBESSELL else if ( BCTAB == 2 ) then do i = 1,NGIRARDI do j = 1,12 DATA(i,j) = GIRARDI(i,j) end do end do NUMDATA = NGIRARDI else if ( BCTAB == 3 ) then do i = 1,NCODE do j = 1,3 DATA(i,j) = CODE(i,j) end do end do NUMDATA = NCODE else print*,"Incorrect BCTAB specified in calling subroutine BCGET." RETURN end if do i = 1,NUMDATA if (DATA(i,CTEFF) < TEFF) LOWTEFF = max(LOWTEFF,DATA(i,CTEFF)) if (DATA(i,CTEFF) >= TEFF) HIGHTEFF= min(HIGHTEFF,DATA(i,CTEFF)) if (DATA(i,CLOGG) < LOGG) LOWLOGG = max(LOWLOGG,DATA(i,CLOGG)) if (DATA(i,CLOGG) >= LOGG) HIGHLOGG= min(HIGHLOGG,DATA(i,CLOGG)) end do do i = 1,NUMDATA if ( DATA(i,CTEFF)==LOWTEFF .and. DATA(i,CLOGG)==LOWLOGG ) then INTDATA(1) = DATA(i,CBC) else if (DATA(i,CTEFF)==LOWTEFF.and.DATA(i,CLOGG)==HIGHLOGG)then INTDATA(2) = DATA(i,CBC) else if (DATA(i,CTEFF)==HIGHTEFF.and.DATA(i,CLOGG)==LOWLOGG)then INTDATA(3) = DATA(i,CBC) else if(DATA(i,CTEFF)==HIGHTEFF.and.DATA(i,CLOGG)==HIGHLOGG)then INTDATA(4) = DATA(i,CBC) end if end do do i = 1,4 if ( INTDATA(i) < -5.0d0 ) then print*,"Could not find BC coefficient ",i print*,"Teff,lowteff,highteff:",Teff,lowteff,highteff print*,"logg,lowlogg,highlogg:",logg,lowlogg,highlogg RETURN end if end do WEIGHTS(1) = (HIGHTEFF-TEFF) / (HIGHTEFF-LOWTEFF) * & (HIGHLOGG-LOGG) / (HIGHLOGG-LOWLOGG) WEIGHTS(2) = (HIGHTEFF-TEFF) / (HIGHTEFF-LOWTEFF) * & (LOGG-LOWLOGG) / (HIGHLOGG-LOWLOGG) WEIGHTS(3) = (TEFF-LOWTEFF) / (HIGHTEFF-LOWTEFF) * & (HIGHLOGG-LOGG) / (HIGHLOGG-LOWLOGG) WEIGHTS(4) = (TEFF-LOWTEFF) / (HIGHTEFF-LOWTEFF) * & (LOGG-LOWLOGG) / (HIGHLOGG-LOWLOGG) BC = 0.0d0 do i = 1,4 BC = BC + WEIGHTS(i)*INTDATA(i) end do END SUBROUTINE BCGET !====================================================================== SUBROUTINE FLOWER (TEFF1,TEFF2,RADIUS1,RADIUS2,LSUN,MBOLSUN, & V,V1,V2,MV,MV1,MV2,DIST1,DIST2,DIST) implicit none real*8 TEFF1,TEFF2,LSUN,MBOLSUN ! IN: Various quanitites real*8 RADIUS1,RADIUS2 ! IN: Radii in solar units real*8 V,V1,V2 ! IN: V mag of system & stars real*8 MV,MV1,MV2 ! OUT: Absolute V magnitudes real*8 DIST,DIST1,DIST2 ! OUT: Distances in parsecs real*8 LOGTEFF,BC1,BC2,BC ! LOCAL: Bolometric correct'ns real*8 PI,SB,RSUN ! LOCAL: Useful constants real*8 MBOL1,MBOL2,LOGL1,LOGL2 ! LOCAL: useful quantities integer i ! LOCAL: Loop counter LOGTEFF = 0.0d0 BC1 = 0.0d0 BC2 = 0.0d0 do i = 1,2 if ( i == 1 ) LOGTEFF = log10(TEFF1) if ( i == 2 ) LOGTEFF = log10(TEFF2) if ( LOGTEFF > 3.9d0 ) then BC = -159648.484d0 + 185570.891d0*LOGTEFF & - 86180.8631d0*LOGTEFF**2 + 19988.9657d0*LOGTEFF**3 & - 2315.54438d0*LOGTEFF**4 + 107.171921d0*LOGTEFF**5 else if ( LOGTEFF < 3.7d0 ) then BC = -18696.2756d0 + 15214.1739d0*LOGTEFF & - 4128.70051d0*LOGTEFF**2 + 373.629476d0*LOGTEFF**3 else BC = 2876.21236d0-3556.15263d0*LOGTEFF+1598.56781d0*LOGTEFF**2 & - 312.299468d0*LOGTEFF**3 + 22.4941524d0*LOGTEFF**4 endif if ( i == 1 ) BC1 = BC if ( i == 2 ) BC2 = BC end do PI = atan(1.0d0) * 4.0d0 ! Pi to the machine precision SB = 5.670400d-8 ! Stefan-Boltzmann constant RSUN = 6.9599d8 ! Radius of Sun from Zombeck DIST1 = -100.0d0 ! Initialise DIST2 = -100.0d0 ! Initialise DIST = -100.0d0 ! Initialise LOGL1 = log10( (4.0d0*PI*SB*(RADIUS1*RSUN)**2*TEFF1**4) / LSUN ) LOGL2 = log10( (4.0d0*PI*SB*(RADIUS2*RSUN)**2*TEFF2**4) / LSUN ) MBOL1 = (-2.5d0 * LOGL1) + MBOLSUN MBOL2 = (-2.5d0 * LOGL2) + MBOLSUN MV1 = MBOL1 - BC1 MV2 = MBOL2 - BC2 MV = -2.5d0 * log10(10.0d0**(-0.4d0*MV1)+10.0d0**(-0.4d0*MV2)) if ( V1 > -99.0d0 ) DIST1 = 10.0d0 ** ( (V1-MV1+5.0d0)/5.0d0 ) if ( V2 > -99.0d0 ) DIST2 = 10.0d0 ** ( (V2-MV2+5.0d0)/5.0d0 ) if ( V > -99.0d0 ) DIST = 10.0d0 ** ( (V-MV+5.0d0)/5.0d0 ) END SUBROUTINE FLOWER !======================================================================