! Program to simulate the Measured Mass absorption efficienies during CARES.
! Written by Pat Arnott. Started 4 August 2013.
! V2 (version 2) of this code chooses aveBC from a specific lognormal distribution function,
! rather than uniformally picking aveBC from a range of values.
! Version 2 created on 9 August 2013.
IMPLICIT NONE
! Declare PAS relevant variables.
REAL aveBabs,BabsTrial,Babs(1000000) ! Perfect Babs, Trial value, and accepted value.
REAL dBabs ! Uncertainty in absorption measurement including accuracy and precision.
REAL dBabsPrecision
REAL BabsAccuracyFraction ! Measurement uncertainty based only on calibration, eg 5%.
REAL LaserPower ! Laser Power used in the photoacoustic instrument, in mW, at the res. freq.
! Declare SP2 relevant variables.
REAL aveBC,BCTrial,BC(1000000) ! Perfect BC, Trial value, and accepted value.
REAL aveBCtrial
REAL BCmin,BCmax ! Range of BC values from the experiment, in ug/m3.
REAL dBC ! Uncertainty in BC mass measurement including accuracy and precision.
REAL dBCPrecision
REAL BCAccuracyFraction ! Measurement uncertainty baed only on calibration, eg 30%.
REAL mu,std ! average value and standard deviation of the log normal distribution function for BC.
REAL mode ! exp(mu - std^2) is the value of the domain where log normal is maximum.
REAL peak ! value of the lognormal distribution evaluated at the mode value.
! Other variables.
REAL TauMax ! Maximum measurement time in seconds for both instruments.
REAL MACperfect ! Slope of the Babs vs BC curve for an ideal experiment.
REAL MAC ! Value of MAC found in a simulation of the experiment.
INTEGER numberOfMeasurements
REAL NormalDis,LogNormalDis,RANDinRANGE ! Define variable type for the functions associated with them.
INTEGER i,j ! i is the loop over all potential values. j counts the accepted values.
REAL P1,P1max,RoleOfDice,fac
! Calculate sqrt(2*pi) as is it used often. Define it as fac.
fac = SQRT(2. * 4.*atan(1.0))
! Set the fixed values.
MACperfect = 21.0 ! m^2 / g
numberOfMeasurements = 48 * 07 ! 60 30 minute averages per day * 14 days of measurements.
! T0 value.
* BCmin = 0.06 ; BCmax = 0.18 ! Range of BC values from the experiment, in ug/m3.
! T1 value.
BCmin = 0.04 ; BCmax = 0.24 ! Range of BC values from the experiment, in ug/m3.
! If you are using T0, comment out the T1 values and remove the comment for the T0 values.
! T0 values for the BC distribution. These are measured values.
! mu = -1.81
! std = 0.41
! T1 values for the BC distribution. These are measured values.
mu = -2.24
std = 0.27
! Compute the mode; needed to choose a value of aveBC.
! These values are fixed by the choice of mu and std.
mode = EXP(mu - std*std)
peak = LogNormalDis(mode,mu,std)
! NOTE: Laser power is different for each photoacoustic instrument.
LaserPower = 20.0 ! mW laser power for the 355 nm instrument.
TauMax = 600.0 ! maximum measurement time in seconds for instrument stable values. 10 minutes.
dBabsPrecision = 8.0 * 30.0 / (LaserPower * SQRT(TauMax))
dBCPrecision = 10.0 / SQRT(TauMax) ! ng/m3 units.
dBCPrecision = dBCPrecision / 1000.0 ! ug/m3 units.
! Instrument calibration uncertainties, related to measurement accuracy.
BabsAccuracyFraction = 0.05 ! PAS calibration uncertainty is +- 5%.
BCAccuracyFraction = 0.30 ! SP2 calibration uncertainty is +- 30%.
! Assume for a given experiment that the instruments are accurate enough.
! BabsAccuracyFraction = 0.00 ! PAS calibration uncertainty is +- 5%.
! BCAccuracyFraction = 0.00 ! SP2 calibration uncertainty is +- 30%.
! Initialize the measurement counter.
j=0
! Create the output file.
OPEN(1,FILE='MACsimulation.txt')
WRITE(1,*) 'MeasNum aveBC_ug/m3 aveBabs_Mm-1 deltaBC_ug/m3 deltaBabs_Mm-1 BC_ug/m3 Babs_Mm-1 MACperfect MAC'
! Test many values to find 'good' values of Babs and BC for making histograms with.
DO 1 i=1,100000000
! Step A). Choose a value for aveBC from the range BCmin to BCmax.
! aveBC is chosen from the measured log normal distribution function for BC.
aveBCtrial=RANDinRANGE(BCmin,BCmax)
P1 = LogNormalDis(aveBCtrial,mu,std)/peak ! Probability normalized to be 1 at its peak value.
IF (RAND() .GT. P1) GOTO 1 ! Reject this trial value of aveBC. Start the process over.
! If the program is here, then a value of aveBC has been accepted. Ready to move to Step B.
aveBC=aveBCtrial
! Step B). Choose a corresponding value of aveBabs.
aveBabs = MACperfect * aveBC ! Choose a corresponding value for Babs.
! Step C). Choose a trial value of Babs.
dBabs = SQRT( dBabsPrecision**2 + (aveBabs*BabsAccuracyFraction)**2 )
BabsTrial = RANDinRANGE(-4.0*dBabs+aveBabs,4.0*dBabs+aveBabs)
! IF (BabsTrial.LT.dBabs) GOTO 1 ! Reject this value of BabsTrial because it is too close to zero, i.e. is deep in the noise.
! Step D). Calculate the probability that the trial value of Babs is good.
P1max = 1.0 / (dBabs * fac)
P1 = NormalDis(BabsTrial,aveBabs,dBabs)
! Step E). Test to see if the trial value of Babs will be accepted or not.
RoleOfDice = RANDinRANGE(0.0,P1max)
IF (RoleOfDice .GT. P1) GOTO 1 ! The trial value of BabsTrial has FAILED.
! Otherwise, if we've made it this far, than the trial value of BabsTrial is good.
! Now generate a trial value of BC and test it.
! Step F) Generate and evaluate a trial value of BC now.
! Calculate the uncertainties in BC.
dBC = SQRT( dBCPrecision**2 + (aveBC * BCAccuracyFraction)**2 )
BCTrial = RANDinRANGE(-4.0*dBC+aveBC,4.0*dBC+aveBC)
! IF (BCTrial.LT.dBC) GOTO 1 ! Reject this value as being within 1 standard deviation of zero, i.e. very uncertain.
P1max = 1.0 / (dBC * fac)
P1 = NormalDis(BCTrial, aveBC, dBC)
RoleOfDice = RANDinRANGE(0.0,P1max)
IF (RoleOfDice .GT. P1) GOTO 1 ! The trial value of BabsTrial has FAILED.
! Step G. We have found good values of BCTrial and BabsTrial. Record them to arrays and output file.
j = j + 1
BC(j) = BCTrial
Babs(j) = BabsTrial
MAC = Babs(j) / BC(j)
WRITE(*,*) 'Trial#, Measurement#, Babs, BC, MAC ',i,j,Babs(j),BC(j),MAC
WRITE(1,*) j,aveBC,aveBabs,dBC,dBabs,BC(j),Babs(j),MACperfect,MAC
! Now end measurements if we have completed the required number.
IF (j.EQ.NumberOfMeasurements) GOTO 2
1 CONTINUE ! Otherwise, go back to 1, and repeat the Do loop again.
2 CLOSE(1)
END
! Function to get a random variable in the range from xmin to xmax.
REAL FUNCTION RANDinRANGE(xmin,xmax)
IMPLICIT NONE
REAL xmin,xmax
RANDinRANGE = xmin + (xmax - xmin)*RAND()
RETURN
END
! Normal Distribution Function.
REAL FUNCTION NormalDis(x,xmean,std)
IMPLICIT NONE
REAL x,xmean,std
NormalDis=0.39894228*exp(-((x-xmean)**2)/(2.0*(std**2)))/std
RETURN
END
! LogNormal Distribution Function.
REAL FUNCTION LogNormalDis(x,xmean,std)
IMPLICIT NONE
REAL x,xmean,std
LogNormalDis=0.39894228*exp(-((LOG(x)-xmean)**2)/(2.0*(std**2)))/(std*x)
RETURN
END