PROGRAM UVBYSPLIT ! John Southworth jkt~astro.keele.ac.uk ! Astrophysics Group, Keele University, UK ! This inputs the overall Stromgren colours and the light ratios of a ! binary star and outputs the individual colours with robust errors. ! Version 1 of the code: 19/11/2004 !======================================================================= implicit none integer NPAR ! number of input parameters integer NRES ! number of output results parameter ( NPAR=8 , NRES=12 ) real*8 PAR(NPAR),PARERR(NPAR) ! input parameters and errors real*8 RES(NRES),RESERR(NRES) ! output results and errors real*8 ERRORS(NPAR,NRES) ! individual error contributions real*8 EHIGH(NRES),ELOW(NRES),SUM ! various stuff integer i,j,ERROR ! useful integers character NAMES(NRES)*10 ! output result names data NAMES / "system V ","system b-y","system m_1","system c_1", & "star:A V ","star:A b-y","star:A m_1","star:A c_1", & "star:B V ","star:B b-y","star:B m_1","star:B c_1"/ 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 write(*,'(A40,A40)') "UVBYSPLIT John Southworth (Keele ", & " University, UK, jkt~astro.keele.ac.uk)" write(*,'(A40,A40)') "This program inputs the Stromgren indic", & "es of a binary star and the light ratios" write(*,'(A40,A40)') "in each passband, and outputs the separ", & "ate Stromgren indices of the two stars," write(*,'(A40,A40)') "with an error budget for each one, expre", & "ssed as percentages of the total errors." write(*,'(A40,A40)') "Enter the Stromgren V b-y m_1 c_1 f", & "or the combined light of the two stars: " read (*,*) PAR(1),PAR(2),PAR(3),PAR(4) write(*,'(A40,A40)') "Enter the absolute errors in these four ", & "quantities: " read (*,*) PARERR(1),PARERR(2),PARERR(3),PARERR(4) write(*,'(A40,A40)') "Enter the light ratios of the two stars ", & "in the u v b y passbands: " read (*,*) PAR(5),PAR(6),PAR(7),PAR(8) write(*,'(A40,A40)') "Enter the absolute errors in these four ", & "quantities: " read (*,*) PARERR(5),PARERR(6),PARERR(7),PARERR(8) ! First calculate the results for the given parameters. Then ! for each parameter, perturb its value by plus and by minus ! its statued 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 ( 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 for each result, add up the individual uncertainties, ! in quadrature, casued by the uncertainty in each parametr. ! The overall uncertainty for each result is the square root ! of the sum of the squares of the individual uncertainties. ! They scale the uncertainties to be expressed as fractions ! of the overall uncertainties for the results calculated. 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 ! Now output the results, their overall uncertainties, units ! and the individual contributions from each parameter. write (*,*) " " write (*,'(A40,A40)') "----------------------------------------", & "------------------------------- " write (*,'(A40,A40)') "Parameter Value Total Fraction", & "al error contributions " write (*,'(A40,A40)') " error V b-", & "y m_1 c_1 lr_u lr_v lr_b lr_y " write (*,'(A40,A40)') "----------------------------------------", & "------------------------------- " do j = 1,NRES if ( RES(j) > -99.0d0 ) then write (*,500) NAMES(j),RES(j),RESERR(j) do i = 1,NPAR 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 (*,*) " " ! End the line end if end do write (*,'(A40,A40)') "----------------------------------------", & "------------------------------- " write (*,*) " " 500 FORMAT (A10," =",F8.4," +/-",F7.4,1X,$) 501 FORMAT (1X,F4.2,$) 502 FORMAT (1X,F4.3,$) 503 FORMAT (A5,$) END PROGRAM UVBYSPLIT !======================================================================= SUBROUTINE CALC (IN,NIN,NOUT,OUT) ! IN(1) is system V ! OUT(3) is system m1 ! IN(2) is system b-y ! OUT(4) is system c1 ! IN(3) is system m1 ! OUT(5) is star 1 V ! IN(4) is system c1 ! OUT(6) is star 1 b-y ! IN(5) is u light ratio ! OUT(7) is star 1 m1 ! IN(6) is v light ratio ! OUT(8) is star 1 c1 ! IN(7) is b light ratio ! OUT(9) is star 2 V ! IN(8) is y light ratio ! OUT(10) is star 2 b-y ! OUT(1) is system V ! OUT(11) is star 2 m1 ! OUT(2) is system b-y ! OUT(12) is star 2 c1 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 U,V,B,Y ! LOCAL: system uvby values real*8 U1,V1,B1,Y1,U2,V2,B2,Y2 ! LOCAL: individual uvby value integer i,j ! LOCAL: loop counters Y = IN(1) B = IN(2) + Y V = IN(3) + 2.0d0*B - Y U = IN(4) + 2.0d0*V - B U1 = U - 2.5d0 * log10( 1.0d0 / (1.0d0+IN(5))) U2 = U - 2.5d0 * log10( IN(5) / (1.0d0+IN(5))) V1 = V - 2.5d0 * log10( 1.0d0 / (1.0d0+IN(6))) V2 = V - 2.5d0 * log10( IN(6) / (1.0d0+IN(6))) B1 = B - 2.5d0 * log10( 1.0d0 / (1.0d0+IN(7))) B2 = B - 2.5d0 * log10( IN(7) / (1.0d0+IN(7))) Y1 = Y - 2.5d0 * log10( 1.0d0 / (1.0d0+IN(8))) Y2 = Y - 2.5d0 * log10( IN(8) / (1.0d0+IN(8))) OUT(1) = IN(1) OUT(2) = IN(2) OUT(3) = IN(3) OUT(4) = IN(4) OUT(5) = Y1 OUT(6) = B1 - Y1 OUT(7) = V1 - 2.0d0*B1 + Y1 OUT(8) = U1 - 2.0d0*V1 + B1 OUT(9) = Y2 OUT(10) = B2 - Y2 OUT(11) = V2 - 2.0d0*B2 + Y2 OUT(12) = U2 - 2.0d0*V2 + B2 END SUBROUTINE CALC !=======================================================================