!======================================================================= ! PROGRAM JKTEBOP John Taylor j.k.taylor~warwick.ac.uk ! Astrophysics Group Warwick University !----------------------------------------------------------------------- ! V(1) = surface brightness ratio V(14) = tidal lead/lag angle (deg) ! V(2) = sum of fractional radii V(15) = third light ! V(3) = ratio of stellar radii V(16) = phase correction ! V(4) = limb darkening of star 1 V(17) = light scaling factor ! V(5) = limb darkening of star 2 V(18) = integration ring size (deg) ! V(6) = orbital inclination V(19) = orbital period (days) ! V(7) = e cos(omega) OR ecentricity V(20) = ephemeris timebase (days) ! V(8) = e sin(omega) OR omega VEXTRA(1) = primary star radius ! V(9) = gravity darkening 1 VEXTRA(2) = secondary star radius ! V(10)= gravity darkening 2 VEXTRA(3) = stellar light ratio ! V(11) = primary reflected light VEXTRA(4) = eccentricity ! V(12) = secondary reflected light VEXTRA(5) = periastron longitude ! V(13) = stellar mass ratio VEXTRA(6) = reduced chi-squared !----------------------------------------------------------------------- ! Version 1: Simplex minimisation algorithm and TASK3 subroutine used ! Version 2: Monte Carlo simulation subroutine added ! Version 3: Adjustments to Monte Carlo LD coeffs and input/output files ! Version 4: Now solves for sum of radii. Convergence criterion modified ! Version 5: Added TASK0 to find LD and GD coeffs. Other minor mods. ! Version 6: Reflections and scale factor can all be fixed or adjusted. ! Version 7: Fit for e,w; SFACT modified; obs'l errors found; spherical. ! Version 8: Modified to perform a bootstrapping error analysis (TASK4). ! Version 9: Command-line args, no-MC-param-kick option, fit for P & T0. ! Version 10: 99999 datapoints, magnitudes used, TASK4+TASK5 => TASK45 ! Last modified: 4th November 2005 (working on psi Centauri). !----------------------------------------------------------------------- ! Possible modifications in future: ! 2) extend to WD2003 and WINK ! 3) allow constraints on light ratio ! 5) incorporate nonlinear limb darkening laws ! 6) remove data 3sigma from best fit, refit, output clipped light curve ! 9) port to F90 or F95 to use long lines, modules, improved output ! 11) JVC suggests that formal fitting errors are x2 too small ! 12) implement parameter kicking to check that minima are global ! 13) incorporate change of omega (apsidal motion) ! 14) add other LD and GD to TASK0 ! 17) go through subroutine LIGHT to understand it all ! 19) add a genetic algorithm for finding robust best fit to light curve !----------------------------------------------------------------------- ! Miscellaneous notes: ! 1) Phase shift is redefined so that primary minimum is moved to the ! phase corresponding to the input value. ! 2) MRQMIN adjusts coeffs only if VARY (called 'ia' in MRQMIN) is 1. ! 3) If VARY=2 then the parameter is fixed during the initial fit but is ! perturbed by a set amount (flat distribution) during Monte Carlo. ! 4) If VARY(11) and/or VARY(12) are -1 then V(11) and/or V(12) are ! calculated from the system geometry; if 0 they are fixed at the ! input value and if 1 are freely adjusted to best fit. ! 5) If the mass ratio is less than or equal to zero then both stars are ! assumed to be spherical. ! 6) If ecosw is greater than 5.0 then the values of ecosw and esinw ! will be taken to be (10+e) and omega and fitting will occur using ! e and omega as parameters. This is a bad idea normally as e and ! omega can be strongly correlated, but can be useful if either e or ! omega are already known and its value is not being adjusted. ! 7) Observational errors are looked for in the input light curve file. ! If they are not found then equal weight is given to each point. !----------------------------------------------------------------------- ! Task numbers and purposes: ! (0) This outputs LD and GD coefficients given Teffs and loggs. ! (1) This outputs a model light curve for fixed input parameters. ! (2) This fits a model to an observed light curve and outputs results. ! (3) This investigates how different parameters vary around best fit. ! (5) This conducts bootstrapping simulations to find robust errors. ! (4) This conducts Monte Carlo simulations to find robust errors. !======================================================================= ! Language: JKTEBOP is written in FORTRAN 77 using standard F77 syntax ! with a few straightforward extras, and compiled using g77. ! Possible non-F77 bits: == <= < > >= /= ! and part-line comments !======================================================================= !======================================================================= PROGRAM JKTEBOP implicit none integer TASK ! Task to perform (between 0 and 5) character INFILE*30 ! Name of the input parameter file real*8 V(20) ! Fundamental EBOP model parameters integer VARY(20) ! Params fixed (0) or variable (1) real*8 DATA(3,99999) ! Data: (time,magnitude,uncertainty) integer NUMDATA ! Number of data points integer ERROR ! Error flag for input/output integer NMONTE ! Number of Monte Carlo to do TASK5 character*1 PERTURB ! whether to perturb simulation pars PERTURB = "y" write (*,*) " " write (*,'(A40,A40)') "JKTEBOP John Taylor (Physics Dept,", & "Warwick Univ, j.k.taylor~warwick.ac.uk)" write (*,'(A40,A40)') "Task 0 outputs limb and gravity darkeni", & "ng coefficients for given Teff and log g" write (*,'(A40,A40)') "Task 1 outputs one model light curve ca", & "lculated using a set of input parameters" write (*,'(A40,A40)') "Task 2 finds the best fit of the model ", & "to observations (internal errors quoted)" write (*,'(A40,A40)') "Task 3 fits observations and finds good", & "ness of fit for several parameter values" write (*,'(A40,A40)') "Task 4 finds robust reliable errors by ", & "analysing with a bootstrapping algorithm" write (*,'(A40,A40)') "Task 5 finds robust errors by analysing", & " with a Monte Carlo simulation algorithm" write (*,'(A40,A40)') "Empty input files can be obtained by ent", & "ering 'newfile', instead of a file name." ! Now check for command-line arguments. If none are present ! then instead prompt for the task to do and the input file. if ( iargc() /= 2 ) then write (*,'(A40,$)') "Enter the task number and file name >> " read (*,*,IOSTAT=ERROR) TASK,INFILE else CALL GETARG (1,INFILE) read (INFILE,*,iostat=ERROR) TASK CALL GETARG (2,INFILE) end if if ( ERROR /= 0 ) then write (*,'(A38,A42)') "## ERROR reading input: it must be one", & " integer followed by a character string. " write (*,*) " " STOP end if if ( TASK < 0 .or. TASK > 5 ) then write (*,'(A38,A42)') "## ERROR: task integer does not corres", & "pond to a task implemeneted in JKTEBOP. " write (*,*) " " ; STOP end if if ( INFILE == "newfile" .or. INFILE == "NEWFILE" ) then CALL NEWFILE (TASK) write (*,*) " " ; STOP end if if ( TASK == 0 ) CALL TASK0 (INFILE) if ( TASK > 0 ) then write (*,*) " " CALL INPUT (TASK,INFILE,V,VARY,DATA,NUMDATA,NMONTE) if ( NMONTE < 0 ) then PERTURB = "n" NMONTE = abs(NMONTE) end if end if if ( TASK == 1 ) CALL TASK1 (V) if ( TASK == 2 ) CALL TASK2 (V,VARY,DATA,NUMDATA) if ( TASK == 3 ) CALL TASK3 (V,VARY,DATA,NUMDATA) if ( TASK == 4 .or. TASK == 5 ) & CALL TASK45 (TASK,V,VARY,DATA,NUMDATA,NMONTE,PERTURB) write (*,*) " " END PROGRAM JKTEBOP !======================================================================= !======================================================================= SUBROUTINE NEWFILE (TASK) ! Makes all empty input files implicit none integer TASK ! Task the file is wanted for character*30 OUTFILE ! Name of output file to make integer NUMCHAR(4) ! Number of character strings integer i,ERROR ! Loop counter and error flag write (*,'(A40,$)') "Enter name of output file to create >> " read(*,*) OUTFILE ; ERROR = 0 CALL OPENFILE (50,"new"," output ",OUTFILE,ERROR) if ( ERROR /= 0 ) STOP write(50,101)"Sum of the radii , Ratio of the radii " write(50,101)"Orbital inclination , Mass ratio of system " write(50,101)"Orbital eccentricity , Periastron longitude " write(50,101)"Gravity darkening star A , Gravity darkening B " write(50,101)"Surface brightness ratio , Third light " write(50,101)"Limb darkening star A , Limb darkening star B " write(50,101)"Reflection effect star A , Reflection for star B " write(50,101)"Phase shift of Min(A) , Light scale factor " if ( TASK == 2 .or. TASK == 3 .or. TASK == 4 ) then write(50,101)"Orbital period of eclipsing binary system (days) " write(50,101)"Reference time of primary minimum (HJD) " end if if ( TASK == 1 ) then write(50,101)"Output file name (continuous character string) " end if if ( TASK == 4 ) then write(50,101)"Number of bootstrapping samples to do " end if if ( TASK == 5 ) then write(50,101)"Number of Monte Carlo simulations to do " end if if ( TASK >= 2 .and. TASK <= 5 ) then write(50,101)"Adjust RADII SUM or RADII RATIO (0 or 1 or 2)" write(50,101)"Adjust INCLINATION or MASSRATIO (0 or 1 or 2)" write(50,101)"Adjust ECCENTRICITY or OMEGA (0 or 1 or 2)" write(50,101)"Adjust GRAVDARK1 or GRAVDARK2 (0 or 1 or 2)" write(50,101)"Adjust SURFACEBRIGHT2 or THIRDLIGHT (0 or 1 or 2)" write(50,101)"Adjust LIMBDARK1 or LIMBDARK2 (0 or 1 or 2)" write(50,101)"Adjust REFLECTION COEFFS 1 and 2 (-1, 0, 1 ,2)" write(50,101)"Adjust PHASESHIFT or SCALE FACTOR (0 or 1 or 2)" write(50,101)"Adjust PERIOD or TZERO (min light) (0 or 1) " write(50,101)"Name of file containing light curve " write(50,101)"Name of output parameter file " end if if ( TASK == 2 ) then write(50,101)"Name of output light curve file " write(50,101)"Name of output model light curve fit file " end if if ( TASK == 4 .or. TASK == 5 ) then write(50,101)"Name of output file of individual fits " end if write (50,*) " " write (50,*) " " write (50,102) "# Enter the appropriate numbers on the l", & "eft-hand side of each line of this file." write (50,*) " " write (50,102) "# Put a negative number for the mass rat", & "io to force the stars to be spherical. " write (50,102) "# The mass ratio will then be irrelevent", & " as it is only used to get deformations." write (50,*) " " write (50,102) "# If e < 10 then e and omega will be ass", & "umed to be e*cos(omega) and e*sin(omega)" write (50,102) "# If e >= 10 then e and omega will be as", & "sumed to be (e+10) and omega (degrees). " write (50,102) "# The first option is in general better ", & "unless eccentricity is larger or fixed. " write (50,*) " " if ( TASK == 1 ) then write (50,102) "# For Task 1 the reflection coefficients", & " are not calculated from the geometry of" write (50,102) "# the system but are taken from the inpu", & "t file. The calculated coefficients are" write (50,102) "# outputted to the terminal so can be pu", & "t in the input file and the task rerun. " write (50,*) " " end if if ( TASK >= 2 .and. TASK <= 5 ) then write (50,102) "# For each adjustable parameter the adju", & "stment integer can be 0 (parameter will" write (50,102) "# be fixed at the input file value), 1 (", & "parameter will be freely adjusted) or 2" write (50,102) "# (Tasks 3, 4, 5; parameter will be fixe", & "d but perturbed during further analysis)" write (50,*) " " write (50,102) "# When fitting a light curve the reflec", & "tion coefficients can be calculated from" write (50,102) "# the system geometry (put -1 for the a", & "djustment integers), held fixed (put 0)" write (50,102) "# or freely adjusted to fit the light cu", & "rve (put 1) - useful for close binaries." write (50,*) " " end if if ( TASK == 4 .or. TASK == 5 ) then write (50,102) "# When doing Monte Carlo or bootstrappin", & "g the starting parameters for each simu-" write (50,102) "# lation are perturbed by a pre-defined ", & "amount to avoid biassing the results. If" write (50,102) "# this is not wanted, put a minus sign i", & "n front of the number of simulations. " write (50,*) " " end if close (50,status="keep") ; write (*,'(A2,A27,I1,A21,A30)') "An", &" empty input file for task ",TASK," has been written to ",OUTFILE 101 FORMAT (29X,A51) 102 FORMAT (A40,A40) END SUBROUTINE NEWFILE !======================================================================= !======================================================================= SUBROUTINE INPUT (TASK,INFILE,V,VARY,DATA,NUMDATA,NMONTE) ! Reads the input file for each task implicit none integer TASK ! IN: Task to read file for character INFILE*30 ! IN: Name of the input file real*8 V(20) ! OUT: Fundamental EBOP params integer VARY(20) ! OUT: Par vary (1), fixed (0) real*8 DATA(3,99999) ! OUT: Data (time, mag, error) integer NUMDATA ! OUT: Number of data points integer NMONTE ! OUT: Number of Monte Carlos character OBSFILE*30 ! LOCAL: Input lightcurve file character LCFILE*30 ! LOCAL: Output lightcurv file character PARAMFILE*30 ! LOCAL: Output parameter file character FITFILE*30 ! LOCAL: Output model fit file integer i,ERROR,STATUS ! LOCAL: Counter & error flags real FMAG,LP,LS ! LOCAL: LIGHT subroutine out real ECOSW,ESINW ! LOCAL: e*cos(w) and e*sin(w) real*8 ARRAY(99999),SELECT ! LOCAL: Help array & FUNCTION ERROR = 0 ; STATUS = 0 CALL OPENFILE (60,"old","input file",INFILE,ERROR) if ( ERROR /= 0 ) STOP do i = 1,20 V(i) = 0.0d0 ; VARY(i) = 0 end do do i = 1,99999 DATA(1,i) = 0.0d0 DATA(2,i) = 0.0d0 DATA(3,i) = 0.0d0 end do !----------------------------------------------------------------------- CALL READFF (60,"RADII SUM ", 0.0d0, 0.8d0, V( 2), & "RADIIRATIO", 0.0d0, 1.0d1, V( 3),STATUS) CALL READFF (60,"INCLNATION",50.0d0, 1.4d2, V( 6), & "MASS RATIO",-1.0d5, 1.0d1, V(13),STATUS) CALL READFF (60,"ECCENTRCTY", 0.0d0,11.0d0, V( 7), & "OMEGA ",-3.6d2, 3.6d2, V( 8),STATUS) CALL READFF (60,"GRAVDARK-1", 0.0d0, 2.0d0, V( 9), & "GRAVDARK-2", 0.0d0, 2.0d0, V(10),STATUS) CALL READFF (60,"SURF-BRT-2", 0.0d0, 2.0d1, V( 1), & "THIRDLIGHT", 0.0d0, 1.0d1, V(15),STATUS) CALL READFF (60,"LIMBDARK-1", 0.0d0, 1.0d0, V( 4), & "LIMBDARK-2", 0.0d0, 1.0d0, V( 5),STATUS) CALL READFF (60,"REFLECTN 1", 0.0d0, 1.0d0, V(11), & "REFLECTN 2", 0.0d0, 1.0d0, V(12),STATUS) CALL READFF (60,"PHASESHIFT",-1.0d0, 1.0d0, V(16), & "SCALEFACTR",-1.0d3, 1.0d3, V(17),STATUS) if ( TASK >= 2 .and. TASK <= 5 ) then CALL READF (60, "PERIOD ", 0.0d0, 1.0d4, V(19),STATUS) CALL READF (60, "TIME-ZERO ", 0.0d0, 3.0d6, V(20),STATUS) end if if ( TASK == 1 ) then CALL READCHAR30 (60,"OUTFILE ",PARAMFILE,STATUS) end if if ( TASK == 4 .or. TASK == 5 ) then read (60,*,iostat=ERROR) NMONTE if ( ERROR /= 0 ) then write (*,'(A36,A44)') "## ERROR reading number of bootstrap", & "ping or Monte Carlo simulations to do. " STATUS = 1 end if end if if ( TASK >= 2 .and. TASK <= 5 ) then CALLREAD2(60,"adj(R1+R2)","adj(R2/R1)",VARY( 2),VARY( 3),STATUS) CALLREAD2(60,"adj(_INC_)","adj(M2/M1)",VARY( 6),VARY(13),STATUS) CALLREAD2(60,"adj(_ECC_)","adj(OMEGA)",VARY( 7),VARY( 8),STATUS) CALLREAD2(60,"adj(_GD1_)","adj(_GD2_)",VARY( 9),VARY(10),STATUS) CALLREAD2(60,"adj(_SB2_)","adj(_L_3_)",VARY( 1),VARY(15),STATUS) CALLREAD2(60,"adj(_LD1_)","adj(_LD2_)",VARY( 4),VARY( 5),STATUS) CALLREAD2(60,"adj(REFL1)","adj(REFL2)",VARY(11),VARY(12),STATUS) CALLREAD2(60,"adj(PSHFT)","adj(SFACT)",VARY(16),VARY(17),STATUS) CALLREAD2(60,"adj(PERIOD","adj(TZERO)",VARY(19),VARY(20),STATUS) CALL READCHAR30 (60,"OBSFILE ",OBSFILE ,STATUS) CALL READCHAR30 (60,"PARAMFILE ",PARAMFILE,STATUS) end if if ( VARY(20) == 1 .and. VARY(16) == 1 ) then write (*,'(A39,A41)') ">> TZERO and PSHIFT cannot both be adju", & "sted so adj(PSHIFT) has been set to zero." VARY(16) = 0 end if if ( TASK == 2 .or. TASK == 4 .or. TASK == 5 ) & CALL READCHAR30 (60,"LCFILE ",LCFILE, STATUS) if ( TASK == 2 ) CALL READCHAR30 (60,"FITFILE ",FITFILE,STATUS) close (unit=60,status="keep") if ( STATUS /= 0 ) GOTO 666 !----------------------------------------------------------------------- CALL OPENFILE (62,"new","parameter ",PARAMFILE,STATUS) write(62,'(A40,A40)') "========================================", & "========================================" write(62,'(A40,A40)') "JKTEBOP output results ", & " John Taylor j.k.taylor~warwick.ac.uk" write(62,'(A40,A40)') "========================================", & "========================================" write(62,*) " " if ( TASK >= 2 .and. TASK <= 5 ) then CALL READDATA (61,OBSFILE,DATA,NUMDATA,STATUS) if ( STATUS /= 0 ) GOTO 666 end if if ( TASK == 2 .or. TASK == 4 .or. TASK == 5 ) & CALL OPENFILE(63,"new","lightcurve",LCFILE, STATUS) if (TASK == 2) CALL OPENFILE(64,"new","model fit ",FITFILE,STATUS) if ( STATUS /= 0 ) GOTO 666 !----------------------------------------------------------------------- if ( V(13) < 0.0d0 ) write (62,'(A19,A61)') "Mass ratio is below", & " zero: the stellar shapes will be forced to be spherical. " if ( V(7) < 5.0d0 ) write (62,'(A19,A61)') "Eccentricity is bel", & "ow 5: [e,omega] are taken to be [e*cos(omega),e*sin(omega)]." if ( V(7) > 5.0d0 ) write (62,'(A19,A61)') "Eccentricity is mor", & "e than 5: e and omega will be taken to be (e+10) and omega. " write (62,*) " " if ( V(7) < 5.0d0 ) then if (V(7)<-1.0d0.or.V(7)>1.0d0.or.V(8)<-1.0d0.or.V(8)>1.0d0) then write (*,'(A36,A44)') "## ERROR: e*sin(omega) or e*cos(omeg", & "a) are unphysical: must be between -1 and 1." STATUS = 1 end if end if ! Find reflection coefficients if they are not fixed CALL LIGHT (V,0.0,FMAG,LP,LS) if (VARY(11) == -1) V(11) = 0.4d0 * LS * (V(2)/(1.0d0+V(3)))**2 if (VARY(12) == -1) V(12)=0.4d0*LP*(V(2)/(1.0d0+(1.0d0/V(3))))**2 if ( TASK >= 2 .and. TASK <= 5 ) then do i = 11,12 100 FORMAT (A4,I2,A37,I1,A36) if ( VARY(i) == -1 ) write (62,100) & "Adj(",i,") = -1 so reflection effect for star ", & i-10," is calculated from system geometry. " if ( VARY(i) == 0 ) write (62,100) & "Adj(",i,") = -1 so reflection effect for star ", & i-10," is fixed at the input file value. " if ( VARY(i) == 1 ) write (62,100) & "Adj(",i,") = -1 so reflection effect for star ", & i-10," is freely adjusted to the best fit. " if ( VARY(i) == 2 ) write (62,100) & "Adj(",i,") = -1 so reflection effect for star ", & i-10," is set to input value but perturbed." end do write (62,*) " " ! Now to avoid problems with partial derivatives later if ( V(7)>5.0d0 .and. VARY(7)+VARY(8)/=0) V(7)=max(V(7),10.001) if ( V(7) < 6.0d0 .and. V(7) >= 1.0d0 ) V(7) = 0.0d0 if ( V(11) < 0.001d0 .and. VARY(11) == 1 ) V(11) = 0.001d0 if ( V(12) < 0.001d0 .and. VARY(12) == 1 ) V(12) = 0.001d0 end if ! Find a good starting value for the light scale factor if ( VARY(17) /= 0 ) then do i = 1,NUMDATA ARRAY(i) = DATA(2,i) end do V(17) = select(ARRAY,NUMDATA,int(0.2d0*NUMDATA)) end if if ( STATUS == 0 ) RETURN 666 write (*,*) " " ! If there were problems then stop now STOP END SUBROUTINE INPUT !======================================================================= !======================================================================= SUBROUTINE READFF (UNIT,NAME1,LOWER1,UPPER1,VALUE1, & NAME2,LOWER2,UPPER2,VALUE2,STATUS) ! This reads in two double precision numbers and checks if ! they are within given bounds. STATUS is set to 1 if error ! occurs and keeps its previous value if no error occurs. implicit none integer UNIT ! IN: unit number to read from character NAME1*10,NAME2*10 ! IN: names of the parameters real*8 LOWER1,UPPER1,LOWER2,UPPER2 ! IN: allowed ranges of values real*8 VALUE1,VALUE2 ! OUT:values of the parameters integer STATUS ! OUT:set to 1 if error occurs integer ERROR ! LOCAL: error flag for input ERROR = 0 read (UNIT,*,IOSTAT=ERROR) VALUE1,VALUE2 if ( ERROR /= 0 ) then write (*,'(A32,A10,A5,A10)') & "## ERROR reading the parameters ",NAME1," and ",NAME2 STATUS = 1 end if if ( VALUE1 < LOWER1 .or. VALUE1 > UPPER1 ) then write (*,'(A23,A10,A4,F12.5)') & "## ERROR: the value of ",NAME1," is ",VALUE1 write (*,'(A21,F12.5,A4,F12.5)') & "The allowed range is ",LOWER1," to ",UPPER1 STATUS = 1 end if if ( VALUE2 < LOWER2 .or. VALUE2 > UPPER2 ) then write (*,'(A23,A10,A4,F12.5)') & "## ERROR: the value of ",NAME2," is ",VALUE2 write (*,'(A21,F12.5,A4,F12.5)') & "The allowed range is ",LOWER2," to ",UPPER2 STATUS = 1 end if END SUBROUTINE READFF !----------------------------------------------------------------------- SUBROUTINE READF (UNIT,NAME,LOWER,UPPER,VALUE,STATUS) ! This reads one double precision number and checks if it is ! within given bounds. STATUS is set to 1 if an error occurs implicit none integer UNIT ! IN: number of unit to read from character NAME*10 ! IN: name of parameter real*8 LOWER,UPPER ! IN: allowed range of values real*8 VALUE ! OUT: value of the parameter integer STATUS ! IN/OUT: set to 1 if error occurs integer ERROR ! LOCAL: error flag for input ERROR=0 read (UNIT,*,IOSTAT=ERROR) VALUE if ( ERROR /= 0 ) then write (*,'(A31,A10)') "## ERROR reading the parameter ",NAME STATUS = 1 end if if ( VALUE < LOWER .or. VALUE > UPPER ) then write (*,'(A23,A10,A4,F12.5)') & "## ERROR: the value of ",NAME," is ",VALUE write (*,'(A21,F12.5,A4,F12.5)') & "The allowed range is ",LOWER," to ",UPPER STATUS = 1 end if END SUBROUTINE READF !----------------------------------------------------------------------- SUBROUTINE READCHAR30 (UNIT,NAME,VALUE,ERROR) ! This reads in a 30-char value. STATUS is set to 1 if an ! error occurs and keeps its previous value if no error. implicit none integer UNIT ! IN: number of unit to read from character NAME*10 ! IN: name of character to be read character VALUE*30 ! OUT: value of the character integer STATUS ! IN/OUT: set to 1 if error occurs integer ERROR ! LOCAL: error flag for input ERROR = 0 read (UNIT,*,IOSTAT=ERROR) VALUE if ( ERROR /= 0 ) then write (*,'(A28,A30)') "## ERROR reading the string ",NAME STATUS = 1 end if END SUBROUTINE READCHAR30 !----------------------------------------------------------------------- SUBROUTINE READ2 (UNIT,NAME1,NAME2,VALUE1,VALUE2,STATUS) ! This reads two integers on one line and tests if both are ! equal to -1, 0, 1 or 2. If not then STATUS is set to 1. ! If there is no error then STATUS keeps its previous value. implicit none integer UNIT ! IN: number of unit to read from character NAME1*10,NAME2*10 ! IN: names of characters to read integer VALUE1,VALUE2 ! OUT: values of the integers integer STATUS ! IN/OUT: set to 1 if error occurs integer ERROR ! LOCAL: error flag for input ERROR = 0 read (UNIT,*,IOSTAT=ERROR) VALUE1,VALUE2 if ( ERROR /= 0 ) then write (*,'(A30,A10,A5,A10)') & "## ERROR reading the integers ",NAME1," and ",NAME2 STATUS = 1 end if if ( VALUE1 < -1 .or. VALUE1 > 2 ) then write(*,'(A32,A10,A4,I6,A28)')"## Parameter adjustment integer " & ,NAME1," is ",VALUE1," and should be -1, 0, 1 or 2" STATUS = 1 end if if ( VALUE2 < -1 .or. VALUE2 > 2 ) then write(*,'(A32,A10,A4,I6,A28)')"## Parameter adjustment integer " & ,NAME2," is ",VALUE2," and should be -1, 0, 1 or 2" STATUS = 1 end if END SUBROUTINE READ2 !======================================================================= SUBROUTINE READDATA (UNIT,OBSFILE,DATA,NUMDATA,STATUS) implicit none integer UNIT ! IN: Unit number to read file from character*30 OBSFILE ! IN: Name of input light curve file real*8 DATA(3,99999) ! OUT: Obs'l data (time, mag, error) integer NUMDATA ! OUT: Number of datapoints read integer STATUS ! IN/OUT: set to 1 if there is error integer i,ERROR,ERRFLAG ! LOCAL: Loop counter + error flags character*200 CHARHELP ! LOCAL: Helper character string real*8 HELP1,HELP2,HELP3 ! LOCAL: Helper variables ERROR = 0 CALL OPENFILE (61,"old","data file ",OBSFILE,ERROR) if ( ERROR /= 0 ) then STATUS = 1 ; RETURN end if read (UNIT,'(A200)',IOSTAT=ERROR) CHARHELP rewind (UNIT) read (CHARHELP,*,iostat=ERROR) HELP1,HELP2,HELP3 if ( ERROR /= 0 ) then read (CHARHELP,*,iostat=ERROR) HELP1,HELP2 if ( ERROR /= 0 ) then write (*,'(A47,A30)') & "## ERROR: cannot understand first line of file ",OBSFILE return end if ERRFLAG = 0 ! There are no observational errors else ERRFLAG = 1 ! There are observational errors end if do i = 1,99999 if ( ERRFLAG == 0 ) then read (UNIT,*,iostat=ERROR) DATA(1,i),DATA(2,i) DATA(3,i) = -1.0d0 else if ( ERRFLAG == 1 ) then read (UNIT,*,iostat=ERROR) DATA(1,i),DATA(2,i),DATA(3,i) end if if ( ERROR /= 0 ) exit end do NUMDATA = i - 1 if ( NUMDATA < 5 ) then write (*,'(A29)') "## ERROR: too few data to fit" STATUS = 1 end if if ( ERRFLAG == 1 ) then write (62,'(A5,I5,A39,A30)') "Read ",NUMDATA, & " datapoints (with errorbars) from file ",OBSFILE write (*,'(A8,I5,A36,A30)') ">> Read ",NUMDATA, & " datapoints (with errors) from file ",OBSFILE else if ( ERRFLAG == 0 ) then write (62,'(A5,I5,A38,A30)') "Read ",NUMDATA, & " datapoints (no error bars) from file ",OBSFILE write (*,'(A8,I5,A37,A30)') ">> Read ",NUMDATA, & " datapoints (no errorbars) from file ",OBSFILE end if write(62,*) " " close (UNIT) END SUBROUTINE READDATA !======================================================================= SUBROUTINE OPENFILE (UNIT,STATE,FILE,FILENAME,STATUS) ! Opens an existing file. STATUS = 1 if the action was not ! successful and left unchanged if the action was successful implicit none integer UNIT ! IN: unit number to open file to character*30 FILENAME ! IN: name of datafile to open character*10 FILE ! IN: identifier of this datafile character*3 STATE ! IN: "old" or "new" integer STATUS ! OUT: indicates success of opening integer ERROR ! LOCAL: error flag for opening file ERROR = 0 OPEN (unit=UNIT,file=FILENAME,status=STATE,iostat=ERROR) if ( ERROR /= 0 ) write (*,'(A17,A3,1X,A10,A6,A30)') & "## ERROR opening ",STATE,FILE," file ",FILENAME if ( ERROR /= 0 ) STATUS = 1 END SUBROUTINE OPENFILE !======================================================================= !======================================================================= SUBROUTINE TASK0 (OUTFILE) ! This task outputs limb darkening ! coefficients for given temperatures and surface gravities. ! Coefficients are read in from tables and bilinear interpo- ! lation is performed. Limb darkening coefficients come from ! 1993AJ....106.2096 van Hamme uvbyUBVRI ! 1995A+AS..110..329 Diaz-Cordoves,Claret,Gimenez uvbyUBV ! 1995A+AS..114..247 Claret,Diaz-Cordoves,Gimenez RIJHK ! 2000A+A...363.1081 Claret uvbyUBVRIJHK implicit none real*8 TEFF(2),LOGG(2) ! Atmospheric parameters character OUTFILE*30 ! file to output to character FILE(7)*50 ! Coefficients file real*8 COEFFS(1000,19) ! Coefficients read from file integer NCOEFF ! Number of datafile lines real*8 BILINEAR ! FUNCTION real*8 COEFF(13) ! Coefficients to output integer VT(7) ! Microturbulent velocity real*8 MOH(7) ! Metallicity [M/H] integer i,j,k,l,m,ERROR ! Loop counters and error flag integer COL(13) data VT / 0,0,2,2,2,0,4 / data MOH / 0.0,0.0,0.0,+0.5,-0.5,0.0,0.0 / data COL / 4,5,6,7,8,9,10,18,19,11,12,13,14/ FILE(1) = "~/home/tabs/LD-vanhamme.dat" FILE(2) = "~/home/tabs/LD-DiazCordoves.dat" FILE(3) = "~/home/tabs/LD-claret2000-vt2mh+0.dat" FILE(4) = "~/home/tabs/LD-claret2000-vt2mh+5.dat" FILE(5) = "~/home/tabs/LD-claret2000-vt2mh-5.dat" FILE(6) = "~/home/tabs/LD-claret2000-vt0mh+0.dat" FILE(7) = "~/home/tabs/LD-claret2000-vt4mh+0.dat" ERROR = 0 CALL OPENFILE (55,"new","outputfile",OUTFILE,ERROR) if ( ERROR /= 0 ) RETURN ! Now the initial quantities have been defined according to ! what is in the files, get the Teff and logg for each star ! from the screen and check that they are within bounds. write (*,'(A40,$)') "Enter the Teff and logg for star A >> " read (*,*,iostat=ERROR) TEFF(1),LOGG(1) if ( ERROR /= 0 ) then write (*,'(A38,A42)') "### ERROR reading input values. They s", & "hould be two floating-point numbers. " RETURN end if write (*,'(A40,$)') "Enter the Teff and logg for star B >> " read (*,*,iostat=ERROR) TEFF(2),LOGG(2) if ( ERROR /= 0 ) then write (*,'(A38,A42)') "### ERROR reading input values. They s", & "hould be two floating-point numbers. " RETURN end if if ( TEFF(1) < 3500.0d0 .or. TEFF(1) > 45000.0d0 ) then write (*,'(A33,A43,F10.3)') "### ERROR: Teff of star A must be", & " between 3500 and 45000 K and is ",TEFF(1) RETURN end if if ( LOGG(1) < 0.0d0 .or. LOGG(1) > 5.0d0 ) then write (*,'(A33,A45,F8.3)') "### ERROR: logg of star A must be", & " between 0.0 and 5.0 (cm/s) and is ",LOGG(1) RETURN end if if ( TEFF(2) < 3500.0d0 .or. TEFF(2) > 45000.0d0 ) then write (*,'(A33,A43,F10.3)') "### ERROR: Teff of star B must be", & " between 3500 and 45000 K and is ",TEFF(2) RETURN end if if ( LOGG(2) < 0.0d0 .or. LOGG(2) > 5.0d0 ) then write (*,'(A33,A45,F8.3)') "### ERROR: logg of star B must be", & " between 0.0 and 5.0 (cm/s) and is ",LOGG(2) RETURN end if ! Now for the big one... Inside, there is a loop over each ! file (over k) which outputs results for one Teff and logg ! and for one of two groups of filter-dependent coefficients ! Outside this is a loop (over j) for the two stars. ! Outside this is a loop (over i) for two different groups ! of filter-dependent limb darkening coefficients. write (55,*) " " write (55,'(A40,A40)') "Below are linear limb darkening coeffici", & "ents found using bilinear interpolation." write (55,'(A40,A40)') "van Hamme (1993AJ....106.2096V) used the", & " Kurucz (1991) ATLAS9 model atmospheres." write (55,'(A40,A40)') "Diaz-Cordoves & Claret (1995A+AS..110..3", & "29D and 1995A+AS..114..247C) ATLAS9 too." write (55,'(A40,A40)') "Claret (2000A+A...363.1081) used ATLAS9 ", & "with several different [M/H] and Vmicro." do i = 1,2 ! Loop over doing different formats do j = 1,2 ! Loop over the two stars if ( i == 1 .and. j == 1 ) then write (55,*) " " write (55,'(A32,A48)') " T_eff log_g Vt M/H ", & " u v b y U B V " else if ( i == 2 .and. j == 1 ) then write (55,*) " " write (55,'(A32,A48)') " T_eff log_g Vt M/H ", & " R_c I_c R_j I_j J H K " end if write (55,'(A32,A48)') "--------------------------------", & "------------------------------------------------" do k = 1,7 ! Loop over each coefficient file open (35,file=FILE(k),status="old",iostat=ERROR) if ( ERROR /= 0 ) write(*,'(A40,A40)') & "## ERROR: could not open existing file ",FILE(k) if ( ERROR /= 0 ) RETURN do l = 1,1000 if ( k == 1 ) read(35,*,iostat=ERROR) (COEFFS(l,m),m=1,19) if ( k >= 2 ) read(35,*,iostat=ERROR) (COEFFS(l,m),m=1,14) if ( ERROR /= 0 ) exit end do NCOEFF = l - 1 do l = 1,13 if ( k == 1 ) then COEFF(l)=BILINEAR(COEFFS,NCOEFF,19,TEFF(j),LOGG(j),COL(l)) else if ( k >= 2 .and. l <= 12 ) then COEFF(l)= BILINEAR(COEFFS,NCOEFF,19,TEFF(j),LOGG(j),l+2) end if end do if ( k == 1 ) then if ( i == 1 ) then write (55,100) int(TEFF(j)),LOGG(j),(COEFF(m),m=1,7) else if ( i == 2) then write (55,101) int(TEFF(j)),LOGG(j),(COEFF(m),m=8,13) end if else if ( k == 2) then if ( i == 1 ) then write (55,102) int(TEFF(j)),LOGG(j),(COEFF(m),m=1,7) else if ( i == 2) then write (55,103) int(TEFF(j)),LOGG(j),(COEFF(m),m=8,12) end if else if ( k >= 3 ) then if ( i == 1 ) write (55,104) & int(TEFF(j)),LOGG(j),VT(k),MOH(k),(COEFF(m),m=1,7) if ( i == 2 ) write (55,105) & int(TEFF(j)),LOGG(j),VT(k),MOH(k),(COEFF(m),m=8,12) end if end do end do end do CLOSE (55) 100 FORMAT ("vanHamme(93) ",I5,1X,F5.3,7X,7(1X,F6.4)) 101 FORMAT ("vanHamme(93) ",I5,1X,F5.3,7X,5(1X,F6.4),8X,F6.4) 102 FORMAT ("DiazCordoves ",I5,1X,F5.3,7X,7(1X,F6.4)) 103 FORMAT ("DiazCordoves ",I5,1X,F5.3,7X,2(1X,F6.4),14X,3(1X,F6.4)) 104 FORMAT ("Claret(2000) ",I5,1X,F5.3,1X,I1,1X,F4.1,7(1X,F6.4)) 105 FORMAT ("Claret(2000) ",I5,1X,F5.3,1X,I1,1X,F4.1, & 2(1X,F6.4),14X,3(1X,F6.4)) END SUBROUTINE TASK0 !======================================================================= DOUBLEPRECISION FUNCTION BILINEAR (ARRAY,NDATA,NCOEFF,IN1,IN2,COL) ! This performs bilinear interpolation on ARRAY to find the ! interpolated value for IN1 and IN2. ARRAY is assumed to ! have the two discrete distributions in the first two ! columns and then NCOEFF further columns, one of which ! (given by COL) contains the results for interpolation. implicit none integer NDATA ! IN: length of array ARRAY integer NCOEFF ! IN: width of array ARRAY real*8 ARRAY(1000,NCOEFF) ! IN: array containing tabular data real*8 IN1,IN2 ! IN: values to interpolate at integer COL ! IN: column containing quantities real*8 LOW1,HIGH1,LOW2,HIGH2 ! LOCAL: adjacent values to IN1, IN2 real*8 COEFFS(4),WT(4) ! LOCAL: coeffs to sum using weights integer i,ERROR ! LOCAL: loop counter and error flag data COEFFS /-1.0d6,-1.0d6,-1.0d6,-1.0d6/ LOW1 = -1.0d6 ; HIGH1 = 1.0d6 LOW2 = -1.0d6 ; HIGH2 = 1.0d6 ! First go through ARRAY to find the values which bracket ! the two values to interpolate to (IN1 and IN2). Avoid ! getting the same values for LOW and HIGH (which would ! occur if input value was exactly the same as a tabulated ! one) by careful choice of inequalities. do i = 1,NDATA if (ARRAY(i,1) < IN1) LOW1 = max(LOW1, ARRAY(i,1)) if (ARRAY(i,1) >= IN1) HIGH1 = min(HIGH1,ARRAY(i,1)) if (ARRAY(i,2) < IN2) LOW2 = max(LOW2, ARRAY(i,2)) if (ARRAY(i,2) >= IN2) HIGH2 = min(HIGH2,ARRAY(i,2)) end do ! Now find the coeffs which are valid for the four values ! which surround the wanted point in 2D space. do i = 1,NDATA if ( ARRAY(i,1) == LOW1 .and. ARRAY(i,2) == LOW2 ) then COEFFS(1) = ARRAY(i,COL) else if ( ARRAY(i,1) == HIGH1 .and. ARRAY(i,2) == LOW2) then COEFFS(2) = ARRAY(i,COL) else if ( ARRAY(i,1) == LOW1 .and. ARRAY(i,2) == HIGH2 ) then COEFFS(3) = ARRAY(i,COL) else if ( ARRAY(i,1) == HIGH1 .and. ARRAY(i,2) == HIGH2) then COEFFS(4) = ARRAY(i,COL) end if end do ! Now calculate the weights for the four values, which ! depend on how close they are to the wanted point. ! Total weights should add up to unity, so check this. WT(1) = (HIGH1-IN1) / (HIGH1-LOW1) * (HIGH2-IN2) / (HIGH2-LOW2) WT(2) = (IN1-LOW1) / (HIGH1-LOW1) * (HIGH2-IN2) / (HIGH2-LOW2) WT(3) = (HIGH1-IN1) / (HIGH1-LOW1) * (IN2-LOW2) / (HIGH2-LOW2) WT(4) = (IN1-LOW1) / (HIGH1-LOW1) * (IN2-LOW2) / (HIGH2-LOW2) BILINEAR = 0.0d0 do i = 1,4 BILINEAR = BILINEAR + WT(i)*COEFFS(i) end do END FUNCTION BILINEAR !======================================================================= SUBROUTINE TASK1 (V) ! Produces a model light curve implicit none real*8 V(20) ! The Fundamental Parameters integer i,j,k,ERROR ! Loop counters and error flag real FMAG,LP,LS ! LIGHT subroutine output real PHASE ! Phase of evaluation CALL LIGHT (V,0.0,FMAG,LP,LS) V(11) = 0.4d0 * LS * (V(2)/(1.0d0+V(3)))**2 V(12) = 0.4d0 * LP * (V(2)/(1.0d0+(1.0d0/V(3))))**2 write (62,'(A20)') "# PHASE MAGNITUDE " do i = 0,10000 PHASE = i / 10000.0d0 CALL LIGHT (V,PHASE,FMAG,LP,LS) write (62,'(F8.6,1X,F12.8)') PHASE,FMAG end do close (UNIT=62) END SUBROUTINE TASK1 !======================================================================= SUBROUTINE TASK2 (V,VARY,DATA,NUMDATA) ! Fit model parameters to an observed light curve implicit none real*8 V(20) ! The Fundamental Parameters integer VARY(20) ! Param vary (1) or fixed (0) real*8 DATA(3,99999) ! Time, magnitude, uncertainty integer NUMDATA ! Number of observations real*8 CHISQ,VERR(20) ! Fit goodness & param errors integer ITER,IFAIL ! Iteration number and success CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,0) CALL FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,IFAIL) CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,1) write (*,'(A32,I3,A45)') ">> Best fit has been found after",ITER, & " iterations and written to the parameter file" END SUBROUTINE TASK2 !======================================================================= SUBROUTINE TASK3 (V,VARY,DATA,NUMDATA) ! Fit model parameters to ! an observed light curve. Then investigate each adjustable ! parameter by fixing it at a range of values round the best ! fit and seeing what the effect is on the other parameters. implicit none real*8 V(20) ! The Fundamental Parameters real*8 DATA(3,99999) ! Time, magnitude, uncertainty integer VARY(20) ! Param vary (1) or fixed (0) integer NUMDATA ! Number of observations integer ITER ! Numberof functionevaluations integer i,j,k,l,ERROR,STATUS,IFAIL ! Loop counters and errorflags integer NUMVARY,VWHERE(20) ! How many / which params vary real*8 VSTORE(20),VTRY ! Store for bestfit parameters real*8 CHISQ,VERR(20) ! Fit goodness & param errors integer VARYFLAG NUMVARY = 0 ; j = 1 do i = 1,20 if ( i < 9 .or. i > 12 ) then if ( VARY(i) /= 0 ) then NUMVARY = NUMVARY + 1 VWHERE(j) = i j = j + 1 end if end if end do CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,0) CALL FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,IFAIL) CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,2) write (*,'(A40,A40)') ">> Fitting process completed. Best fit i", & "s found and outputted to parameter file. " do i = 1,20 VSTORE(i) = V(i) end do VARYFLAG = 0 do i = 1,NUMVARY do j = 1,20 V(j) = VSTORE(j) end do j = VWHERE(i) if ( j /= 16 .and. j /= 17 ) then if ( VARY(i) == 1 ) VARYFLAG = 1 if ( VARY(i) == 2 ) VARYFLAG = 2 VARY(j) = 0 CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,100+j) if ( j == 6 ) then do k = 1,21 V(6) = VSTORE(6) - 1.0d0 + 0.1d0*(k-1) CALL FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,IFAIL) CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,200+j) end do else V(j) = VSTORE(j) * 0.9d0 CALL FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,IFAIL) do k = 1,21 V(j) = VSTORE(j) * (0.8 + 0.02d0*(k-1)) CALL FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,IFAIL) CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,200+j) end do end if if ( VARYFLAG == 1 ) VARY(i) = 1 if ( VARYFLAG == 2 ) VARY(i) = 2 end if end do END SUBROUTINE TASK3 !======================================================================= SUBROUTINE TASK45 (TASK,V,VARY,DATA,NUMDATA,NMONTE,PERTURB) ! also does the bootstrapping now ! Fits model parameters to an observed light curve and esti- ! mates realistic parameter uncertainties using Monte Carlo ! simulations. This is an excellent error indicator but ! explicitly assumes that there is negligible difference ! between the derived light variation and the actual light ! variation of the binary. NMONTE simulations occur where ! the model light curve is resampled at the phases of the ! observations, Gaussian errors (of ERRSIZE times the actual ! rms ofthe residuals of the best fit to the observed light ! curve) are added and the synthetic light curve is refitted. ! This method can give optimistic results when the best fit ! is not a good fit; bootstrapping must then be used instead. ! The final output contains 1 sigma uncertainties. implicit none integer TASK ! 4 bootstrapping, 5 Monte Carlo real*8 V(20) ! The fundamental parameters real*8 VEXTRA(6) ! Extra (dependent) parameters real*8 DATA(3,99999) ! Time, magnitude and uncertainty integer NMONTE ! Number of simulations to do real*8 VERR(20) ! Parameter errors from fit function integer VARY(20),NVARY ! Param vary (1) or fixed (0) integer NUMDATA ! Number of observations integer ITER ! Number of function iteration integer i,j,k,ERROR,STATUS ! Loop counters and error flags real*8 VALL(26,NMONTE) ! All the best-fitting params found integer SEEDSTART,SEED ! FUNCTION start random number maker real*8 RANDOMG ! FUNCTION for Gaussian random noise real*8 RANDOM ! FUNCTION flat-distrib random noise real*8 DV(20) ! Perturbations to apply to paramtrs real*8 ARRAY(NMONTE) ! All results for one parameter real*8 HELP1,HELP2,HELP3 ! Useful variable storage real FMAG,LP,LS ! LIGHT subroutine variables real PHASE,A,B,EPSILON ! Calculated dependent parameters real EPSILON1,EPSILON2 ! Calculated dependent parameters real*8 INDATA(3,99999) ! Synthetic light curve to be fitted real*8 INV(20) ! Parameters of synthetic lightcurve real*8 CHISQ,ERRSIZE ! Chi-squared of the best fit real*8 PHASER,SELECT,SIGMA ! FUNCTIONS character PERTURB*1 ! Define parameter labels for output and the amounts to vary ! each by if they have not been adjusted to the best fit but ! are required to be perturbed during Monte Carlo analysis. DATA DV / 0.05d0,0.1d0, 0.1d0,-0.1d0,-0.1d0, 0.0d0, 0.1d0, & 0.1d0, 0.1d0, 0.1d0, 0.1d0, 0.1d0, 0.1d0, 0.1d0, & 0.5d0, 0.1d0, 0.1d0, 0.0d0, 1.d-6,-1.d-4 / ! Now find the best fit to the light curve and output it. ! Store some dependent quantities for output later on. CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,0) CALL FITEBOP(DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,ERROR) CALL OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,2) write (*,'(A40,A40)') ">> Fitting process completed. Best fit i", & "s found and outputted to parameter file. " CALL LIGHT (V,0.25,FMAG,LP,LS) VEXTRA(1) = V(2) / (1.0d0 + V(3)) ! r_1 VEXTRA(2) = V(2) / (1.0d0 + (1.0d0/V(3))) ! r_2 VEXTRA(3) = dble(LS / LP) ! L_2 / L_1 VEXTRA(4) = sqrt(V(7)**2 + V(8)**2) ! e VEXTRA(5) = atan2(V(8),V(7))*57.2957795d0 ! omega if (VEXTRA(5)<0.0d0) VEXTRA(5)=VEXTRA(5)+360.0d0 ! check omega VEXTRA(6) = CHISQ write (63,'(A)') "# N ITER SB2 r_1+r_2 "// & " k LD1 LD2 i "// & " e(cos)w e(sin)w GD1 GD2 "// & " refl1 refl2 q tidalangle "// & " L_3 phase_corr sfact int_ring "// & " P T_zero r_1 r_2 "// & " L_2 / L_1 e omega (chi)^2 " NVARY = 0 ! Calculate the number of variable params do i = 1,20 if ( VARY(i) == 1 ) NVARY = NVARY + 1 end do do i = 1,NUMDATA INDATA(1,i) = DATA(1,i) INDATA(2,i) = DATA(2,i) INDATA(3,i) = DATA(3,i) end do ! If the input NMONTE is below zero no perturbations should ! be applied to the initial light curve parameter estimates ! prior to each MC simulation. This is useful if convergence ! only happens for a very narrow range of parameter values. if (PERTURB == "n") write (62,'(A20,A60)') "Number of simulation", & "s is below zero: params will not be perturbed before each. " if (PERTURB == "y") write (62,'(A20,A60)') "Number of simulation", & "s is above zero: params will be perturbed before each one. " write (62,*) " " ! For Monte Carlo simulations must evaluate the best-fitting ! model light curve at the observed phases. Also, if there ! were no errors in the input light curve, must calculate ! the rms of the residuals of the fit and use this to scale ! the Gaussian noise added to the simulated light curve. if ( TASK == 5 ) then HELP1 = 0.0d0 do i = 1,NUMDATA PHASE = real(phaser(DATA(1,i),V(19),V(20))) CALL LIGHT (V,PHASE,FMAG,LP,LS) HELP1 = HELP1 + (DATA(2,i) - FMAG)**2 DATA(2,i) = dble( FMAG ) end do ERRSIZE = sqrt(HELP1 / (NUMDATA-NVARY) ) if ( DATA(3,1) <= 0.0d0 ) then do i = 1,NUMDATA INDATA(3,i) = ERRSIZE end do write (62,'(A36,A44)') "No observational errors were supplie", & "d with the input light curve. Noise of size" write (62,'(F6.3,A14,A60)')abs(ERRSIZE)*1.d3,"mmag (standard", & " error of best fit) will be added to the synthetic datasets." else write (62,'(A36,A44)') "Observational errors were supplied w", & "ith the input light curve. These have been " write (62,'(A36,A44)') "assumed to be correct and used to se", & "t the size of the simulated Gaussian noise. " end if end if SEED = SEEDSTART () !Start random number generator do i = 1,20 INV(i) = V(i) ! Preserve best-fit parameter end do ! values and copy to INV. ! Now add Gaussian noise to the theoretical curve sampled at ! the observed phases. Then perturb the parameters using a ! flat distribution of noise (and put them in INV). Fit the ! synthetic data, reject non-convergent fits, ensure that ! inclination is below 90 degrees and output the results. if ( TASK == 4 ) write(*,'(A45,I6)') & ">> Number of bootstrapping simulations to do:",NMONTE if ( TASK == 5 ) write(*,'(A43,I6)') & ">> Number of Monte Carlo simulations to do:",NMONTE write(*,'(A13,$)') ">> Completed:" !----------------------------------------------------------------------- do i = 1,NMONTE ! Look over simulations 500 continue ! For bootstrapping, randomly sample the observations (with ! replacement) to create a new light curve to fit. ! For Monte Carlo, create a simulated light curve by adding ! Gaussian noise to the best-fit model light curve. if ( TASK == 4 ) then do j = 1,NUMDATA k = 1 + int( random(SEED) * (real(NUMDATA)-0.000001) ) INDATA(1,j) = DATA(1,k) INDATA(2,j) = DATA(2,k) INDATA(3,j) = DATA(3,k) end do else if ( TASK == 5 ) then do j = 1,NUMDATA INDATA(2,j) = DATA(2,j) + randomg(SEED,0.0d0,INDATA(3,j)) end do end if ! Now perturb the initial values of the fitted parameters do j = 1,20 if ( (VARY(j)==1 .and. PERTURB=="y") .or. VARY(j) == 2 ) then if ( DV(j) > 0.0d0 ) then INV(j) = V(j) - 0.5d0*V(j)*DV(j) + RANDOM(SEED)*V(j)*DV(j) else if ( DV(j) <= 0.0d0 ) then INV(j) = V(j) - 0.5d0*DV(j) + RANDOM(SEED)*DV(j) end if end if end do if ( INV(6) > 90.0d0 ) GOTO 500 ! Avoid inclination > 90 ! Now fit the simulated light curve and store the results CALL FITEBOP (INDATA,NUMDATA,INV,VARY,ITER,CHISQ,VERR,ERROR) if ( ERROR /= 0 ) GOTO 500 do j = 1,20 VALL(j,i) = INV(j) end do CALL LIGHT (INV,0.25,FMAG,LP,LS) VALL(21,i) = INV(2) / (1.0d0 + INV(3)) ! r_1 VALL(22,i) = INV(2) / (1.0d0 + (1.0d0/INV(3))) ! r_2 VALL(23,i) = dble(LS / LP) ! L_2 / L_1 VALL(24,i) = sqrt(INV(7)**2 + INV(8)**2) ! e VALL(25,i) = atan2(INV(8),INV(7))*57.2957795d0 ! omega if(VALL(25,i)<0.0d0) VALL(25,i)=VALL(25,i)+360.0d0 !check omega VALL(26,i) = CHISQ ! chisq write(63,'(I5,1X,I3,1X,30(f14.6,1X))') & i,ITER,(INV(j),j=1,20),(VALL(j,i),j=21,26) if ( i < 10 ) write (*,'(1X,I1,$)') i if ( (real(i)/10.0) == int(real(i)/10.0) ) then if ( i >= 10 .and. i < 100 ) write (*,'(1X,I2,$)') i if ( i >= 100 .and. i < 1000 ) write (*,'(1X,I3,$)') i if ( i >= 1000 .and. i < 10000 ) write (*,'(1X,I4,$)') i if ( i >= 10000 ) write (*,'(1X,I5,$)') i end if CALL STOPCHECK (ERROR) if ( ERROR /= 0 ) then NMONTE = i exit end if end do write (*,*) " " CALL OUTPUTSIM (TASK,V,VARY,VEXTRA,VALL,NMONTE) END SUBROUTINE TASK45 !======================================================================= SUBROUTINE STOPCHECK (STATUS) ! Checks if the file "jktebop.stop" ! exists and if so puts STATUS=1000 implicit none integer ERROR,STATUS STATUS = 0 ERROR = 100 open (99,file="jktebop.stop",status="old",iostat=ERROR) close (99,status="delete") if ( ERROR == 0 ) STATUS = 1000 END SUBROUTINE STOPCHECK !======================================================================= SUBROUTINE OUTPUTSIM (TASK,V,VARY,VEXTRA,VALL,NMONTE) implicit none real*8 V(20) ! The fundamental parameters real*8 VEXTRA(6) ! Extra (dependent) parameters integer NMONTE ! Number of simulations to do integer i,j,k,ERROR,STATUS ! Loop counters and error flags real*8 VALL(26,NMONTE) ! All the best-fitting params found real*8 ARRAY(NMONTE) ! All results for one parameter real*8 HELP1,HELP2,HELP3 ! Useful variable storage character*5 NAME(26) ! Names of the variables, quantities real*8 PHASER,SELECT,SIGMA ! FUNCTIONS real*8 PERCENT integer VARY(20),TASK DATA NAME/" J","r1+r2"," k"," u_1"," u_2"," inc","ecosw", & "esinw"," GD1"," GD2","refl1","refl2"," q","T_lag", & " L_3"," d(P)","sfact","iring","P_orb"," T_0 "," r_1", & " r_2","L2/L1"," e","omega","Rchi2"/ ! Now for each adjusted parameter (and the other calculated ! quantities, put the NMONTE parameter evaluations into an ! array and select the median and 1_sigma bounds (using the ! median +/- 0.5_sigma) and output to the parameter file. write (62,*) " " write (62,'(A40,A40)') "========================================", & "========================================" if ( TASK == 4 ) write (62,'(A13,I6,A27)') & "Results from ",abs(NMONTE)," bootstrapping simulations:" if ( TASK == 5 ) write (62,'(A13,I6,A25)') & "Results from ",NMONTE," Monte Carlo simulations:" write (62,'(A40,A40)') "========================================", & "========================================" write (62,*) " " write (62,'(A40,A40)') " Best fit stand.dev. m", & "edian +68.3% -68.3% % " do j = 1,26 do k = 1,NMONTE ARRAY(k) = VALL(j,k) end do HELP1=SELECT(ARRAY,NMONTE,int(NMONTE*0.1585d0)) ! Get median and HELP2=SELECT(ARRAY,NMONTE,int(NMONTE*0.5d0)) ! 68.3% (1sigma) HELP3=SELECT(ARRAY,NMONTE,int(NMONTE*0.8415d0)) ! interval. PERCENT = 50.d0 * abs((abs(HELP3-HELP2)+abs(HELP1-HELP2))/HELP2) if ( j < 20 .and. VARY(j) /= 0 ) then if ( j == 1 .or. j == 2 .or. j == 3 ) then write (62,100) j,NAME(j),mod(V(j),1.0d3),SIGMA(ARRAY,NMONTE) & ,mod(HELP2,1.0d3),HELP2-HELP1,HELP3-HELP2,PERCENT else if ( j /= 14 .and. j /= 17 ) then write (62,101) j,NAME(j),mod(V(j),1.0d3),SIGMA(ARRAY,NMONTE) & ,mod(HELP2,1.0d3),HELP2-HELP1,HELP3-HELP2 end if end if if ( j == 21 .or. j == 22 .or. j == 23 ) then write (62,102) NAME(j),VEXTRA(j-20),SIGMA(ARRAY,NMONTE), & HELP2,HELP2-HELP1,HELP3-HELP2,PERCENT else if ( j == 24 .or. j == 25 .or. j == 26) then write (62,103) NAME(j),VEXTRA(j-20),SIGMA(ARRAY,NMONTE), & HELP2,HELP2-HELP1,HELP3-HELP2 end if end do write (62,*) " " 100 FORMAT (I2,1X,A5,1X,F11.7," +/-",F10.7,2X, & F11.7," +",F10.7," -",F10.7," (",F5.2,"%)") 101 FORMAT (I2,1X,A5,1X,F11.7," +/-",F10.7,2X, & F11.7," +",F10.7," -",F10.7) 102 FORMAT (3X,A5,1X,F11.7," +/-",F10.7,2X, & F11.7," +",F10.7," -",F10.7," (",F5.2,"%)") 103 FORMAT(3X,A5,1X,F11.7," +/-",F10.7,2X,F11.7," +",F10.7," -",F10.7) END SUBROUTINE OUTPUTSIM !======================================================================= SUBROUTINE OUTPUT (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,WHAT) ! This subroutine performs the output for TASK 2 and it is ! called by subroutine TASK2. It has three aspects: output ! of the starting position (L<0), output of the results of ! each iteration (L>0) and final output on convergence (L=0) implicit none real*8 DATA(3,99999),V(20),VERR(20) ! Observations, params & error integer NUMDATA,NVARY ! Number of obs & varying pars integer VARY(20),VWHERE(20) ! Params vary & their position integer ITER ! Num. of function evaluations integer WHAT ! Signifies what to output character VNAME(20)*18,VNAM(24)*5 ! Parameter names, short names integer i,j,k,ERROR,STATUS ! Loop counters and errorflags real PHASE ! Output quantities for data real A,B,EPSILON1,EPSILON2,EPSILON ! Axes, & stellar oblatenesses real FMAG,LP,LS ! LIGHT subroutine variables real*8 RESIDSUM,RESIDABSSUM,CHISQ ! Summing variables for the ..reduced chisq real*8 RESIDSQSUM,FLUXRESIDSQSUM ! .. magnitude, flux residuals real*8 ECC,OMEGA,ECOSW,ESINW ! Orbital shape parameters real*8 PHASER ! FUNCTION to phase datapoints real*8 RP,RS ! Fractional radii ofthe stars real*8 LTOTAL ! Total system light (inc. L3) real*8 SE ! Standard error of one obsv'n DATA VNAME/ "Surface br ratio ","Sum of the radii ", & "Ratio of the radii","Limb darkening 1 ","Limb darkening 2 ", & "Orbit inclination ","ecc cos(omega) ","ecc sin(omega) ", & "Gravity dark 1 ","Gravity dark 2 ","Reflected light 1 ", & "Reflected light 2 ","Mass ratio q ","Tidal lead/lag ang", & "Third light ","Phase correction ","Light scale factor", & "Integration ring ","Orbital period ","Ephemeris timebase"/ DATA VNAM/" J ","r1+r2"," k "," u_1 "," u_2 "," inc ","ecosw", & "esinw"," GD1 "," GD2 ","refl1","refl2"," q ","T_lag", & " L_3 "," d(P)","sfact","iring","P_orb"," T_0 "," r_1 ", & " r_2 ","L2/L1"," rms "/ ! If WHAT=0 output initial parameters and useful quantities ! If WHAT=1 then output the final fitted parameters, useful ! quantities, residual curve, and best-fitting model curve ! If WHAT=2 output best fit parameters and useful quantities ! If WHAT=101-116 then this has been invoked from TASK3 to ! output the name of the parameter (WHAT-100) ! If WHAT=201-216 then this has come from TASK3 to output ! fit results for a certain value of parameter (WHAT-200) RP = V(2) / (1.0d0 + V(3)) ! Find primary radius RS = V(2) / (1.0d0 + (1.0d0/V(3))) ! Find secondary radius CALL LIGHT (V,0.0,FMAG,LP,LS) if ( VARY(11) == -1 ) V(11) = 0.4d0 * LS * RP**2 if ( VARY(12) == -1 ) V(12) = 0.4d0 * LP * RS**2 LTOTAL = LP + LS + V(15) NVARY = 0 do i = 1,20 if ( VARY(i) == 1 ) NVARY = NVARY + 1 end do !----------------------------------------------------------------------- ! First output the current values of the parameters ! Also include information on the dependent variables ! (such as oblateness, secondary radius and light ratio). if ( WHAT == 0 ) then write(62,107) "Initial values of parameters: " else if ( WHAT == 1 .or. WHAT == 2 ) then write(62,*) " " write(62,107) "--------------------------------------- " write(62,110) ITER," iterations of EBOP completed " write(62,107) "--------------------------------------- " write(62,*) " " write(62,107) "Final values of parameters: " end if if ( WHAT == 0 .or. WHAT == 1 .or. WHAT == 2 ) then if ( V(6) > 90.0d0 ) V(6) = 180.0d0 - V(6) if ( V(6) < 0.0d0 ) V(6) = V(6) + 180.0d0 if ( WHAT == 0 ) then do i = 1,20 if ( i == 7 .and. V(7) > 5.0d0 ) then if (VARY(i)/=1) write(62,100) i,"Eccentricity ", & V(i)-10.0d0 if (VARY(i)==1) write(62,101) i,"Eccentricity ", & V(i)-10.0d0," (adjusted)" else if ( i == 8 .and. V(7) > 5.0d0 ) then if (VARY(i)/=1) write(62,100) i,"Periast. longitude",V(i) if (VARY(i)==1) write(62,101) i,"Periast. longitude",V(i), & " (adjusted)" else if ( i /= 18 ) then if (VARY(i)/=1) write(62,100) i,VNAME(i),V(i) if (VARY(i)==1)write(62,101)i,VNAME(i),V(i)," (adjusted)" end if end do else do i = 1,20 if ( i == 7 .and. V(7) > 5.0d0 ) then if ( VARY(i) /= 1 ) write (62,113) i,"Eccentricity ", & V(i)-10.0d0," (fixed)" if ( VARY(i) == 1 ) write (62,114) i,"Eccentricity ", & V(i)-10.0d0,VERR(i) else if ( i == 8 .and. V(7) > 5.0d0 ) then if ( VARY(i) /= 1 ) write (62,113) i,"Periast. longitude", & V(i)-10.0d0," (fixed)" if ( VARY(i) == 1 ) write (62,114) i,"Periast. longitude", & V(i)-10.0d0,VERR(i) else if ( i /= 18 ) then if (VARY(i)/=1) write(62,113) i,VNAME(i),V(i)," (fixed)" if (VARY(i)==1) write(62,114) i,VNAME(i),V(i),VERR(i) end if end do end if write (62,*) " " write(62,102)"Fractional primary radius: ",RP write(62,102)"Fractional secondary radius: ",RS if ( V(2) >= 0.6d0 ) then write (62,*) " " write (62,105) "## WARNING: average radius is greater th", & "an 0.3 so the radii may be wrong by 5%!!" write (62,105) "## See North & Zahn (2004, New.Ast.Rev.,", & " 47, 741) for details. Do not use EBOP!!" write (62,*) " " else if ( V(2) >= 0.5d0 ) then write (62,*) " " write (62,105) "## WARNING: average radius is greater th", & "an 0.25 so the radii may be wrong by 1%." write (62,105) "## See North & Zahn (2004, New Astronomy", & " Review, 47, 741) for further details. " write (62,*) " " end if write(62,102)"Stellar light ratio (at phase 0.25):",LS/LP write(62,102)"Primary contribut'n to system light:",LP/LTOTAL write(62,102)"Secondary contrib'n to system light:",LS/LTOTAL write(62,*) " " if ( V(7) > 5.0d0 ) then ECOSW = (V(7)-10.0d0) * cos(V(8)/57.2957795d0) ESINW = (V(7)-10.0d0) * sin(V(8)/57.2957795d0) write (62,102) "Eccentricity*cos(omega): ",ECOSW write (62,102) "Eccentricity*sin(omega): ",ESINW else if ( V(7) /= 0.0d0 .and. V(8) /= 0.0d0 ) then ECC = sqrt(V(7)**2 + V(8)**2) OMEGA = atan2(V(8),V(7)) * 45.0d0 / atan(1.0d0) if ( OMEGA < 0.0d0 ) OMEGA = OMEGA + 360.0d0 write (62,102) "Eccentricity: ",ECC write (62,102) "Omega (degrees): ",OMEGA end if CALL BIAX(real(RP),real(abs(V(13))),A,B,EPSILON1) CALL BIAX(real(RS),real(1.0/abs(V(13))),A,B,EPSILON2) if ( V(13) <= 0.0d0 ) EPSILON1 = 0.0d0 if ( V(13) <= 0.0d0 ) EPSILON2 = 0.0d0 write (62,102) "Oblateness of the primary star: ",EPSILON1 write (62,102) "Oblateness of the secondary star: ",EPSILON2 if ( (WHAT == 1 .or. WHAT == 2) .and. V(13) <= 0.0d0 ) then CALL BIAX(real(RP),real(abs(V(13))),A,B,EPSILON) write(62,102)"Expected oblateness of primary: ",EPSILON CALL BIAX(real(RS),real(1./abs(V(13))),A,B,EPSILON) write(62,102)"Expected oblateness of secondary: ",EPSILON end if if ( EPSILON1 > 0.04 .or. EPSILON2 > 0.04 ) then write (62,*) " " write (62,105) "## WARNING: oblateness is above the reco", & "mmended maximum value for EBOP of 0.04. " write (62,105) "See Popper & Etzel (1981, AJ, 86, 102) f", & "or a justification and further details. " end if write (62,*) " " end if !----------------------------------------------------------------------- ! Now output the observations and the calculated residuals. if ( WHAT == 1 .or. WHAT == 2 ) then if ( WHAT == 1 ) write (63,'(A62)') & "# TIME MAGNITUDE ERROR PHASE MODEL (O-C) " RESIDSUM = 0.0d0 RESIDABSSUM = 0.0d0 RESIDSQSUM = 0.0d0 do i = 1,NUMDATA PHASE = real( phaser(DATA(1,i),V(19),V(20)) ) CALL LIGHT (V,PHASE,FMAG,LP,LS) if ( WHAT == 1 ) write (63,104) DATA(1,i),DATA(2,i),DATA(3,i), & PHASE,FMAG,(DATA(2,i)-FMAG) RESIDSUM = RESIDSUM + DATA(2,i) - FMAG RESIDABSSUM = RESIDABSSUM + abs(DATA(2,i)-FMAG) RESIDSQSUM = RESIDSQSUM + (DATA(2,i)-FMAG)**2 end do write(62,102) "Sum of the residuals (magnitudes): ", RESIDSUM write(62,102) "Sum of absolute values of residuals:",RESIDABSSUM write(62,102) "Sum of the squares of the residuals:", RESIDSQSUM end if write(62,115) "Number datapoints / degs of freedom:",NUMDATA,NVARY if ( WHAT == 1 .or. WHAT == 2 ) then write (62,102) "Standard error (one obs'n, mmag): ", & sqrt(RESIDSQSUM / (NUMDATA-NVARY)) * 1000.0d0 write (62,102) "rms of residuals (one obs'n, mmag): ", & sqrt(RESIDSQSUM / NUMDATA) * 1000.0d0 end if write (62,102) "Reduced chi-squared from errorbars: ", CHISQ write (62,*) if ( WHAT == 1 ) then write (64,'(A22)') "# PHASE MAGNITUDE " do i = 0,10000 PHASE = i / 10000.0d0 CALL LIGHT (V,PHASE,FMAG,LP,LS) write (64,'(F8.6,1X,F12.8)') PHASE,FMAG end do end if !----------------------------------------------------------------------- if ( WHAT > 100 .and. WHAT < 117 ) then j = 0 do i = 1,20 if ( VARY(i) == 1 ) then j = j + 1 VWHERE(j) = i end if end do write (62,*) " " write (62,'(A41,A20)') & "Now investigating the fitting parameter: ",VNAME(WHAT-100) write (62,'(A4,2X,24(A5,6X))') "ITER",VNAM(WHAT-100),VNAM(24), & (VNAM(VWHERE(k)),k=1,j),VNAM(21),VNAM(22),VNAM(23) write(*,'(A44,A20)') & ">> Now investigating the fitting parameter: ",VNAME(WHAT-100) end if !----------------------------------------------------------------------- if ( WHAT > 200 .and. WHAT < 217 ) then RESIDSQSUM = 0.0d0 do i = 1,NUMDATA PHASE = real( phaser(DATA(1,i),V(19),V(20)) ) CALL LIGHT (V,PHASE,FMAG,LP,LS) RESIDSQSUM = RESIDSQSUM + (DATA(2,i)-FMAG)**2 end do SE = sqrt(RESIDSQSUM / NUMDATA) * 1000.0d0 j = 0 do i = 1,20 if ( VARY(i) == 1 ) then j = j + 1 VWHERE(j) = i end if end do CALL LIGHT (V,0.25,FMAG,LP,LS) write(62,'(I3,23(F11.7))') ITER,V(WHAT-200),SE,(V(VWHERE(k)), & k=1,j),V(2)/(1.0d0+V(3)),V(2)/(1.0d0+(1.0d0/V(3))),LS/LP end if !----------------------------------------------------------------------- 100 FORMAT (I2,2X,A18,2X,F16.8) 101 FORMAT (I2,2X,A18,2X,F16.8,A12) 102 FORMAT (A36,2X,F14.8) 1021 FORMAT (A36,2X,E14.7) 104 FORMAT (F14.6,1X,F9.4,1X,F8.4,3X,F9.7,1X,F8.4,1X,F7.4) 105 FORMAT (A40,A40) 107 FORMAT (A48) 108 FORMAT (F14.6,F10.6) 110 FORMAT (I6,A34) 112 FORMAT (A39,F7.4) 113 FORMAT (I2,2X,A18,2X,F16.8,A9) 114 FORMAT (I2,2X,A18,2X,F16.8," +/- ",F12.8) 115 FORMAT (A36,7X,I5,1X,I3) END SUBROUTINE OUTPUT !======================================================================= !======================================================================= DOUBLEPRECISION FUNCTION PHASER (HJD,PERIOD,TZERO) implicit none real*8 HJD,PERIOD,TZERO PHASER = (HJD - TZERO) / PERIOD PHASER = PHASER - int(PHASER) if ( PHASER < 0.0d0 ) PHASER = PHASER + 1.0d0 END FUNCTION PHASER !======================================================================= !======================================================================= SUBROUTINE FITEBOP (DATA,NUMDATA,V,VARY,ITER,CHISQ,VERR,ERROR) ! This subroutine uses the MRQMIN algorithm (Press, 1992) to ! find the best-fitting EBOP model (with errors) to the data implicit none real*8 DATA(3,99999) ! IN: Data (time, mag, error) integer NUMDATA ! IN: Number of datapoints real*8 V(20) ! IN: EBOP parameter values integer VARY(20) ! IN: Whether parameters vary integer ITER ! OUT: Number of iterations real*8 CHISQ ! OUT: Reduced chi-squared real*8 VERR(20) ! OUT: Formal parameter errors integer ERROR ! OUT: Whether fit successful real*8 X(99999),Y(99999),SIG(99999) ! LOCAL: data to go intoMRQMIN real*8 LAMBDA ! LOCAL: Marquardt lambda real*8 COVAR(20,20),ALPHA(20,20) ! LOCAL: Covariance etc matrix real*8 OCHISQ ! LOCAL: Previous chi-squared real FMAG,LP,LS ! LOCAL: comes from EBOP code integer i,j,NVARY ! LOCAL: helpful integers CALL LIGHT (V,0.0,FMAG,LP,LS) if (VARY(11) == -1) V(11) = 0.4d0 * LS * (V(2)/(1.0d0+V(3)))**2 if (VARY(12) == -1) V(12)=0.4d0*LP*(V(2)/(1.0d0+(1.0d0/V(3))))**2 do i = 1,NUMDATA X(i) = DATA(1,i) Y(i) = DATA(2,i) SIG(i) = abs(DATA(3,i)) end do LAMBDA = -1.0d0 CALL MRQMIN (X,Y,SIG,NUMDATA,V,VARY,20,COVAR,ALPHA,20,CHISQ, & LAMBDA,ERROR) OCHISQ = 1.0d10 do i = 1,200 if ( V(6) > 89.9d0 .and. VARY(6) /= 0) V(6) = 89.9d0 CALL MRQMIN (X,Y,SIG,NUMDATA,V,VARY,20,COVAR,ALPHA,20,CHISQ, & LAMBDA,ERROR) if ( LAMBDA >= 1.0d7 .and. abs(OCHISQ/CHISQ) < 1.001 ) exit if ( ERROR /= 0 ) exit OCHISQ = CHISQ end do ITER = i if ( ERROR /= 0 ) write (60,'(20(F15.6,1X))') V if ( ERROR /= 0 ) RETURN LAMBDA = 0.0d0 CALL MRQMIN (X,Y,SIG,NUMDATA,V,VARY,20,COVAR,ALPHA,20,CHISQ, & LAMBDA,ERROR) do i = 1,20 VERR(i) = sqrt(COVAR(i,i)) end do if ( V(6) > 90.0d0 ) V(6) = 180.0d0 - V(6) NVARY = 0 do i = 1,20 if ( VARY(i) == 1 ) NVARY = NVARY + 1 end do if ( DATA(3,1) > 0.0d0 ) CHISQ = CHISQ / (NUMDATA-NVARY) if ( DATA(3,1) <= 0.0d0 ) CHISQ = 1.0d0 END SUBROUTINE FITEBOP !======================================================================= SUBROUTINE EBOP (X,V,Y,DYDA,NUMCOEFFS,VARY) ! This evaluates the goodness of fit and the numerical deri- ! vatives for each coefficient for one datapoint using EBOP. implicit none real*8 X,Y ! IN+OUT: Input time, output flux real*8 V(20) ! IN: Parameters for the EBOP fit integer NUMCOEFFS ! IN: Total number of parameters integer VARY(20) ! IN: Which params are being adjusted real*8 DYDA(20) ! OUT: Numerical differentials real PHASE ! LOCAL: Phase of the datapoint X integer i,j,k,ERROR ! LOCAL: helpful integers real FMAG,LP,LS ! LOCAL: EBOP subroutine output real*8 STORE ! LOCAL: to help find derivatives real FMAG1,FMAG2 ! LOCAL: to help find derivatives real*8 DV(20),DVV ! LOCAL: amount to perturb params by real*8 phaser data DV/ 0.01d0, 0.01d0, 0.01d0,-0.05d0,-0.05d0,-0.1d0,-0.001d0, & -0.001d0, 0.1d0, 0.1d0, 0.01d0, 0.01d0, 0.01d0,-1.0d0, & -0.01d0,-0.0001d0,-0.01d0, 0.5d0, 0.1d-6,-0.001d0 / PHASE = real(phaser(X,V(19),V(20))) CALL LIGHT (V,PHASE,FMAG,LP,LS) Y = dble(FMAG) if (VARY(11)==-1) V(11) = 0.4d0*LS*(V(2)/(1.0d0+V(3)))**2 if (VARY(12)==-1) V(12) = 0.4d0*LP*(V(2)/(1.0d0+(1.0d0/V(3))))**2 do i = 1,NUMCOEFFS if ( VARY(i) == 1) then if ( DV(i) >= 0.0d0 ) then DVV = abs(DV(i) * V(i)) else DVV = abs(DV(i)) end if STORE = V(i) V(i) = STORE + DVV PHASE = real(phaser(X,V(19),V(20))) CALL LIGHT (V,PHASE,FMAG1,LP,LS) V(i) = STORE - DVV PHASE = real(phaser(X,V(19),V(20))) CALL LIGHT (V,PHASE,FMAG2,LP,LS) V(i) = STORE DYDA(i) = (FMAG1 - FMAG2) / (2.0d0 * DVV) else DYDA(i) = 0.0d0 end if end do if ( V(6) > 89.8d0) then DYDA(6) = DYDA(6) + sign( (exp(V(6)-90.0d0))**20 , DYDA(6)) end if END SUBROUTINE EBOP !======================================================================= !======================================================================= !======================================================================= SUBROUTINE BIAX (R,Q,A,B,EPS) ! EBOP subroutine to calculate biaxial ellipsoid dimensions ! and oblateness for each star after Chandrasekhar (1933). if ( Q <= 0.0 ) then A = R B = R EPS = 0.0 else A = R * ( 1.0 + (1.0 + 7.0*Q)/6.0 * R**3.0) B = R * ( 1.0 + (1.0 - 2.0*Q)/6.0 * R**3.0) EPS = (A - B) / A B=( (1.0 - EPS) * R**3.0) ** (1.0/3.0) A = B / (1.0 - EPS) end if END SUBROUTINE BIAX !======================================================================= SUBROUTINE GETEW (ECOSW,ESINW,E,W) ! EBOP subroutine to calculate e and w from e(cos)w e(sin)w if ( ECOSW == 0.0 .and. ESINW == 0.0 ) then E = 0.0 W = 0.0 else W = atan2( ESINW,ECOSW ) E = sqrt( ESINW*ESINW + ECOSW*ECOSW ) W = W * 180.0 / 3.141592654 end if END SUBROUTINE GETEW !======================================================================= SUBROUTINE LIGHT (V,PHASE,FMAG,LP,LS) C real*8 V(20) REAL LP,LS,LECL,LE C COMMON /VARIB/ BS,RP,RATIO,UP,US,FI,ECOSW,ESINW,YP,YS,SP,SS,Q, C $ TANGL,EL,DPH,SFACT,DGAM,EXTRA(13) C REF. THE ASTROPHYSICAL JOURNAL, 174 617-628,1972 JUNE 15 C ECLIPSING-BINARY SOLUTIONS BY SEQUENTIAL OPTIMIZATION C OF THE PARAMETERS C C DATA PI,TWOPI,RAD /3.14159265E0,6.28318531E0,0.0174532925E0/ C C DETERMINE PRIMARY AND SECONDARY BIAXIAL DIMENSIONS C USE SPERICAL RADII FOR THE COMPUTATION OF ECLIPSE FUNCTIONS C USE OBLATENESSES FOR THE COMPUTATION OF THE OUTSIDE ECLIPSE C PHOTOMETRIC VARIATIONS WITH LIMB AND GRAVITY DARKENING C BS = real(V( 1)) RP = real(V(2)/(1.0d0+V(3))) RATIO = real(V( 3)) UP = real(V( 4)) US = real(V( 5)) FI = real(V( 6)) if ( V(7) > 5.0d0 ) then ECOSW = real( (V(7)-10.0d0) * cos(V(8)/57.2957795d0) ) ESINW = real( (V(7)-10.0d0) * sin(V(8)/57.2957795d0) ) else ECOSW = real(V( 7)) ESINW = real(V( 8)) end if YP = real(V( 9)) YS = real(V(10)) SP = real(V(11)) SS = real(V(12)) Q = real(V(13)) TANGL = real(V(14)) EL = real(V(15)) DPH = 1.0 - real(V(16)) SFACT = real(V(17))!10.0d0**(-0.4*V(17))) DGAM = 5.0 if ( Q <= 0.0 ) then CALL BIAX (RP,0.0,RPA,RPB,EP) RS=RP*RATIO CALL BIAX (RS,0.0,RSA,RSB,ES) else CALL BIAX (RP,Q,RPA,RPB,EP) RS=RP*RATIO CALL BIAX (RS,1./Q,RSA,RSB,ES) end if C C CORRECT THE OBSERVED PHASE FOR ANY EPOCH ERROR IN EPHEMERIS C THETA=PHASE+DPH C SINI = SIN(FI*RAD) SINI2 = SINI*SINI COSI2 = 1.0E0 - SINI2 C C TRANSLATE TIDAL LEAD/LAG ANGLE TO RADIANS TANGR=TANGL*RAD C C EQUATION 9 C CONVERT PHASE TO RADIANS FMN=THETA*TWOPI C C GET CURRENT VALUES OF E, AND W CALL GETEW (ECOSW,ESINW,E,W) C C TEST FOR CIRCULAR ORBIT IF (E) 17,20,17 20 COSVW=COS(FMN) SINVW=SIN(FMN) RV=1.0E0 GO TO 25 C C SOLUTION OF KEPLER'S EQUATION BY DIFFERENTIAL CORRECTIONS C (NON-ZERO ECCENTRICITY ONLY . . . ) C C EQUATION 6 C 17 OMEGA = 450.0E0 - W 23 IF (OMEGA - 360.0E0) 22,21,21 21 OMEGA = OMEGA - 360.0E0 GO TO 23 22 OMEGA = OMEGA*RAD C SINE AND COSINE OF OMEGA COSW=COS(OMEGA) SINW=SIN(OMEGA) C C COMPUTE MEAN ANOMALY CORRECTION TO PHASE C CORRESPONDING TO V=OMEGA=90-W C AT WHICH PHASE COS(V-OMEGA)=1 E0=ATAN2(SQRT(1.0E0-E*E)*SINW,COSW+E) C C MEAN ANOMALY OF MID-PRIMARY ECLIPSE FMA0=E0-E*SIN(E0) C C MEAN ANOMALY FMA=FMN+FMA0 C FIRST APPROXIMATION OF ECCENTRIC ANOMALY EA=FMA+E*SIN(FMA) C DO 10 J=1,15 C EVALUATE SINE AND COSINE OF ECCENTRIC ANOMALY SINE=SIN(EA) COSE=COS(EA) DENOM=1.0E0-E*COSE DISC=FMA-EA+E*SINE EA=EA+DISC/DENOM C TEST FOR CONVERGENCE IF (ABS(DISC) - 2.0E-05) 15,15,10 10 CONTINUE C C C EVALUATE SINE AND COSINE OF TRUE ANOMALY 15 COSV=(COSE-E)/DENOM SINV=SINE*SQRT(1.0E0-E*E)/DENOM C C RADIUS VECTOR RV = (1.0E0-E*E)/(1.0E0+E*COSV) C C THE PHOTOMETRIC PHASE ARGUMENT IN TERMS OF ORBIT PARAMETERS C VW = V-OMEGA COSVW=COSV*COSW+SINV*SINW SINVW=SINV*COSW-COSV*SINW C 25 COS2=COSVW*COSVW SIN2=1.0E0-COS2 C CSVWT=COS(TANGR)*COSVW-SIN(TANGR)*SINVW C C C PHOTOMETRIC EFFECTS C C C TEST FOR SIMPLE CASE OF TWO SPHERICAL STARS IF (EP .EQ. 0. .AND. ES .EQ. 0.) GO TO 26 C C EITHER OR BOTH STARS ARE OBLATE C FMAXP=((1.0E0-UP)+0.666666667E0*UP*(1.0E0+0.2E0*EP)) 1 *(1.0E0+3.0E0*YP*EP)/(1.0E0-EP) FMAXS=((1.0E0-US)+0.666666667E0*US*(1.0E0+0.2E0*ES)) 1 *(1.0E0+3.0E0*YS*ES)/(1.0E0-ES) C CHANGE IN INTENSITY RATIO DUE TO OBLATENESS RELATED VARIABLES C FROM QUADRATURE TO MINIMUM C FACE ON TO END ON DELTP=(15.0E0+UP)/(15.0E0-5.0E0*UP)*(1.0E0+YP)*EP DELTS=(15.0E0+US)/(15.0E0-5.0E0*US)*(1.0E0+YS)*ES C FORE-SHORTENING FUNCTION OF OBLATENESS SHORT=SINI2*CSVWT*CSVWT GO TO 27 C C BOTH STARS ARE SPHERICAL C 26 FMAXP=1.0E0-UP/3.0E0 FMAXS=1.0E0-US/3.0E0 DELTP=0.0E0 DELTS=0.0E0 SHORT=0.0 C C UN-NORMALIZED BRIGHTNESS OF STELLAR COMPONENTS AT QUADRATURE 27 OP=PI*RPB*RPB*FMAXP OS=PI*RSB*RSB*FMAXS*BS C THE NORMALIZING FACTOR OTOT=OP+OS C BRIGHTNESS CONTRIBUTION FROM EACH COMPONENT LP=OP/OTOT*(1.0E0-DELTP*SHORT) LS=OS/OTOT*(1.0E0-DELTS*SHORT) C C REFLECTION AND RERADIATION EQUATION IF (SP .EQ. 0.0E0 .AND. SS .EQ. 0.0E0) GO TO 28 HEAT=SINI*COSVW HEAT2=0.5E0+0.5E0*HEAT*HEAT DLP=SP*(HEAT2+HEAT) DLS=SS*(HEAT2-HEAT) GO TO 29 28 DLP=0.0E0 DLS=0.0E0 C C WHICH ECLIPSE COULD THIS BE 29 IF (COSVW) 40,40,30 C C PRIMARY ECLIPSE C 30 R1 = RP R2 = RS UU = UP LE=LP DLE=DLP GO TO 60 C C C SECONDARY ECLIPSE C 40 R1 = RS R2 = RP UU = US LE=LS DLE=DLS C 60 SUM = 0.0E0 ALAST = 0.0E0 AREA=0.0E0 C C EQUATION 5 C DD = SINVW*SINVW + COSVW*COSVW*COSI2 IF (DD .LE. 1.0E-06) DD=0.0 DD = DD*RV*RV D = SQRT(ABS(DD)) R22 = R2*R2 C C EQUATION 17 C GAMN = 90.01E0*RAD DGAMA = DGAM*RAD DGM = DGAMA/2.0E0 RK = 0.0E0 GAM = 0.0E0 50 GAM = GAM + DGAMA C HAS LIMIT OF INTEGRATION BEEN REACHED IF (GAM - GAMN) 48,48,49 C 48 RR = R1*SIN(GAM) R12 = RR*RR C AA = 0.0E0 C ARE THE PROJECTED DISKS CONCENTRIC IF (D) 405,406,405 406 IF (RR - R2) 230,230,403 403 IF (RK - R2) 404, 49, 49 404 AA = PI*R22 GO TO 215 C TEST FOR NO ECLIPSE 405 IF (D-R1-R2) 240,216,216 216 SUM = 0.0E0 GO TO 49 C DECIDE WHICH AREA EQUATIONS FOR NON-CONCENTRIC ECLIPSE 240 IF (D-RR-R2) 245,215,215 245 IF (D-R2+RR) 230,230,250 250 IF (R1-R2) 255,255,280 255 IF (DD-R22+R12) 205,210,210 280 IF (D-RR+R2) 290,260,260 260 IF (RR-R2) 255,255,265 265 IF (DD-R12+R22) 270,210,210 C C EQUATION 12 C 270 S1 = ABS((R12 - R22 - DD)*0.5E0/D) A1 = ABS(R2-S1) B2 = ABS(RR-S1-D ) AA=PI*R22-(R22*ACOS((R2-A1)/R2) 1 - (R2-A1)*SQRT(2.0E0*R2*A1-A1*A1)) 2 +R12*ACOS((RR-B2)/RR)-(RR-B2)*SQRT(2.0E0*RR*B2-B2*B2) GO TO 215 C 290 IF (R1 - R2 - D) 260,260,295 295 IF (RK - R2 - D) 300,215,215 300 RR = R2 + D R12 = RR*RR GAMN = 0.0E0 GO TO 260 C 230 AA = PI*R12 GO TO 215 C C EQUATION 10 C 205 S = ABS((R12 - R22 + DD)*0.5E0/D) A = ABS(RR-S) B1 = ABS(R2-S-D) A1 = R12*ACOS((RR-A)/RR) - (RR-A)*SQRT(2.0E0*RR*A - A*A) AB1 = R22*ACOS((R2-B1)/R2) - (R2-B1)*SQRT(2.0E0*R2*B1-B1*B1) AA = PI*R12 - A1 + AB1 GO TO 215 C C EQUATION 1 C 210 S = ABS((R12 - R22 + DD)*0.5E0/D) A = ABS(RR-S) B = ABS(S-D+R2) A1 = R12*ACOS((RR-A)/RR) - (RR-A)*SQRT(2.0E0*RR*A - A*A) AA1 = R22*ACOS((R2-B)/R2) - (R2-B)*SQRT(2.0E0*R2*B - B*B) AA = A1 + AA1 C 215 DAREA = AA - ALAST SUM = SUM + DAREA*(1.0E0 - UU + UU*COS(GAM-DGM)) ALAST = AA AREA = AREA + DAREA C RK = RR GO TO 50 C C LIGHT LOSS FROM ECLIPSE C 49 ADISK = PI*R1*R1 ALPHA = SUM/(ADISK*(1.0E0-UU/3.0E0)) LECL = ALPHA*LE AREA = AREA/ADISK REFL=DLP+DLS-AREA*DLE C C THEORETICAL INTENSITY WITH THIRD LIGHT AND QUADRATURE C SCALE FACTOR APPLIED C C FLITE = ((LP+LS-LECL+REFL)*(1.0E0-EL)+EL)*SFACT FLITE = ((LP+LS-LECL+REFL)*(1.0E0-EL)+EL) FMAG = -2.5d0 * log10(FLITE) + SFACT C RETURN END !======================================================================= !======================================================================= !================= NUMERICAL RECIPES SUBROUTINES =============== !======================================================================= !======================================================================= SUBROUTINE MRQMIN (x,y,sig,ndata,a,ia,ma,covar,alpha,nca,chisq, & alamda,ifail) implicit none integer NDATA,MA,NCA ! IN: numdata, numcoeffs, maxcoeffs real*8 X(ndata),Y(ndata) ! IN: data to be fitted real*8 SIG(ndata) ! IN: data errors in y real*8 A(ma) ! IN: coefficients integer IA(ma) ! IN: adjust (1) or fix (0) coeffs real*8 COVAR(nca,nca) ! OUT: curvature matrix real*8 ALPHA(nca,nca) ! OUT: covariance matrix real*8 ALAMDA ! IN/OUT: Marquardt lambda factor real*8 CHISQ ! OUT: chi-squared of the fit integer MMAX,j,k,l,m,mfit,ifail parameter (MMAX = 20) real*8 ochisq,atry(MMAX),beta(MMAX),da(MMAX) SAVE ochisq,atry,beta,da,mfit if(alamda.lt.0.0d0)then mfit=0 do 11 j=1,ma if (ia(j).eq.1) mfit=mfit+1 11 continue alamda=0.0010d0 call mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nca,chisq) ochisq=chisq do 12 j=1,ma atry(j)=a(j) 12 continue endif j=0 do 14 l=1,ma if(ia(l).eq.1) then j=j+1 k=0 do 13 m=1,ma if(ia(m).eq.1) then k=k+1 covar(j,k)=alpha(j,k) endif 13 continue covar(j,j)=alpha(j,j)*(1.0d0+alamda) da(j)=beta(j) endif 14 continue call gaussj(covar,mfit,nca,da,1,1,ifail) if ( ifail /= 0 ) return if(alamda.eq.0.0d0)then call covsrt(covar,nca,ma,ia,mfit) return endif j=0 do 15 l=1,ma if(ia(l).eq.1) then j=j+1 atry(l)=a(l)+da(j) endif 15 continue call mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,nca,chisq) if(chisq.lt.ochisq)then alamda=0.1d0*alamda ochisq=chisq j=0 do 17 l=1,ma if(ia(l).eq.1) then j=j+1 k=0 do 16 m=1,ma if(ia(m).eq.1) then k=k+1 alpha(j,k)=covar(j,k) endif 16 continue beta(j)=da(j) a(l)=atry(l) endif 17 continue else alamda=10.0d0*alamda chisq=ochisq endif return END !----------------------------------------------------------------------- SUBROUTINE mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nalp,chisq) IMPLICIT NONE integer ma,nalp,ndata,ia(ma),MMAX,mfit,i,j,k,l,m parameter ( MMAX = 20 ) real*8 chisq,a(ma),alpha(nalp,nalp),beta(ma),sig(ndata),x(ndata) real*8 y(ndata),dy,sig2i,wt,ymod,dyda(MMAX) mfit=0 do 11 j=1,ma if (ia(j).eq.1) mfit=mfit+1 11 continue do 13 j=1,mfit do 12 k=1,j alpha(j,k)=0.0d0 12 continue beta(j)=0.0d0 13 continue c print*,alpha;print*," ";print*,beta chisq=0. do 16 i=1,ndata CALL EBOP (x(i),a,ymod,dyda,ma,ia) sig2i=1./(sig(i)*sig(i)) dy=y(i)-ymod j=0 do 15 l=1,ma if(ia(l).eq.1) then j=j+1 wt=dyda(l)*sig2i k=0 do 14 m=1,l if(ia(m).eq.1) then k=k+1 alpha(j,k)=alpha(j,k)+wt*dyda(m) endif 14 continue beta(j)=beta(j)+dy*wt endif 15 continue chisq=chisq+dy*dy*sig2i 16 continue do 18 j=2,mfit do 17 k=1,j-1 alpha(k,j)=alpha(j,k) 17 continue 18 continue return END !----------------------------------------------------------------------- SUBROUTINE GAUSSJ (a,n,np,b,m,mp,ifail) implicit none integer m,mp,n,np,NMAX,ifail real*8 a(np,np),b(np,mp),big,dum,pivinv parameter (NMAX=20) integer i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX) ifail=0 ; irow=0 ; icol=0 do j=1,n ipiv(j)=0 end do do i=1,n big=0.0d0 do j=1,n if(ipiv(j).ne.1)then do k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big)then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then print*,"Singular matrix in gaussj ",k ifail = 1 ; return endif end do endif end do ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum end do do l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum end do endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.0d0) then print*,"singular matrix in gaussj ",ICOL ifail = 1 ; return end if pivinv=1.0d0/a(icol,icol) a(icol,icol)=1.0d0 do l=1,n a(icol,l)=a(icol,l)*pivinv end do do l=1,m b(icol,l)=b(icol,l)*pivinv end do do ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0.0d0 do l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum end do do l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum end do endif end do end do do l=n,1,-1 if(indxr(l).ne.indxc(l))then do k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum end do endif end do return END !----------------------------------------------------------------------- SUBROUTINE COVSRT (covar,npc,ma,ia,mfit) implicit none integer ma,mfit,npc,ia(ma),i,j,k real*8 covar(npc,npc),swap do 12 i=mfit+1,ma do 11 j=1,i covar(i,j)=0.0d0 covar(j,i)=0.0d0 11 continue 12 continue k=mfit do 15 j=ma,1,-1 if(ia(j).eq.1)then do 13 i=1,ma swap=covar(i,k) covar(i,k)=covar(i,j) covar(i,j)=swap 13 continue do 14 i=1,ma swap=covar(k,i) covar(k,i)=covar(j,i) covar(j,i)=swap 14 continue k=k-1 endif 15 continue return END !======================================================================= !======================================================================= !====================================================================== INTEGER FUNCTION SEEDSTART () ! This uses the Fortran-intrinsic function SYSTEM_CLOCK to ! generate an integer SEED dependent on real time. ! SYSTEM_CLOCK returns the present time, the number of ! times a second this changes (100 on my computer) and the ! maximum value of present time (after which present time ! starts again from zero) (2147483647 on my computer). ! SEED is outputted as a four-figure integer. implicit none integer a,b ! unused constant values CALL SYSTEM_CLOCK (SEEDSTART,a,b) ! SEEDSTART becomes the current clock count END FUNCTION SEEDSTART !===================================================================== DOUBLEPRECISION FUNCTION RANDOM (SEED) ! Minimal random number generator from Park and Miller. ! Return a random deviate between 0.0 and 1.0 with a flat ! ditribution. Its period is (2**31 - 1) = 2.15e9 ! The constants are the best found by Park and Miller. ! Set and reset SEED to any integer value to inititalize ! the sequence and do not change it thereafter - the code ! updates SEED itself. ! Does not work if SEED=0, unless it bit-swaps. implicit none integer SEED ! All-important seed value integer IA,IM,IR,IQ,MASK ! Various constants real*8 AM ! Another constant integer K ! Yet another constant IA = 16807 IM = 2147483647 ! 2**31 ie largest integer the .. AM = 1.0d0 / IM ! .. computer can handle. IQ = 127773 IR = 2836 MASK = 123459876 SEED = IEOR(SEED,MASK) ! IEOR is integer exclusive-or K = SEED / IQ ! .. bit-swapping, the non- SEED = IA * (SEED - K*IQ) - K*IR ! .. numeric operation if (SEED < 0.0d0) SEED = SEED + IM ! .. needed by all random RANDOM = AM * SEED ! .. number generstors. SEED = IEOR(SEED,MASK) ! Seed is updated here END FUNCTION RANDOM !======================================================================= DOUBLEPRECISION FUNCTION RANDOMG (SEED,MEAN,SD) ! This produces a random number with a Gaussian distribut- ! ion. It uses RANDG to generate a random number distrib ! with zero mean and unity variance then scales the result ! with SD and offsets it with MEAN to produce the effect. implicit none integer SEED ! Seeding integer real*8 MEAN,SD ! Desired mean and S.D. (variance) integer ISET ! LOCAL: integer variable real*8 FAC,GSET,RSQ,V1,V2 ! LOCAL: real variables real*8 RANDOM ! FUNCTION: flat-distrib random number generator SAVE ISET,GSET if ( SEED < 0 ) ISET = 0 ! Reinitialise if ( ISET == 0 ) then do V1 = 2.0d0 * RANDOM(SEED) - 1.0d0 V2 = 2.0d0 * RANDOM(SEED) - 1.0d0 RSQ = V1**2 + V2**2 if ( RSQ /= 0.0d0 .and. RSQ <= 1.0d0 ) exit end do FAC = SQRT( -2.0d0 * log(RSQ) / RSQ ) GSET = V1 * FAC RANDOMG = V2 * FAC ISET = 1 else RANDOMG = GSET ISET = 0 end if RANDOMG = ( RANDOMG * SD ) + MEAN END FUNCTION RANDOMG !===================================================================== DOUBLEPRECISION FUNCTION SELECT (ARRAY,NUM,K) ! Returns the Kth smallest value in ARRAY(1:NUM). ! ARRAY is simply sorted during this procedure. implicit none integer NUM ! IN: size of the array real*8 ARRAY(NUM) ! IN: Input array integer K ! OUT: value to find real*8 STORE ! LOCAL: storage variable integer i,j,TAG ! LOCAL: loop counters and a flag TAG = 0 do i = 1,NUM STORE = 10.0**10.0 do j = i,NUM if ( ARRAY(j) < STORE ) then STORE = ARRAY(j) TAG = j end if end do ARRAY(TAG) = ARRAY(i) ARRAY(i) = STORE end do SELECT = ARRAY(K) END FUNCTION SELECT !======================================================================= DOUBLEPRECISION FUNCTION SIGMA (ARRAY,NUM) ! Returns the standard deviation of array ARRAY implicit none integer NUM ! IN: size of the array real*8 ARRAY(NUM) ! IN: array real*8 MEAN,VAR,SUM,SUMSQ ! LOCAL: properties of ARRAY values integer i ! LOCAL: loop counter SIGMA = 0.0d0 SUM = 0.0d0 SUMSQ = 0.0d0 do i = 1,NUM SUM = SUM + ARRAY(i) SUMSQ = SUMSQ + ARRAY(i)**2 end do MEAN = SUM / NUM VAR = (SUMSQ / NUM) - MEAN**2 if ( VAR > 0.0d0 ) SIGMA = sqrt( VAR ) END FUNCTION SIGMA !======================================================================= !======================================================================= !======================================================================= !=======================================================================