FORTRAN source code for Fake Data.

This routine is designed to simulated the output of the APDs at the COAST telescope under a variety of experimental conditions. The expectation value of intensity follows a Gaussian modulated sinusoid, as would be expected from a Gaussian shaped optical filter. Photon shot noise is simulated by selecting the photon counts from a Poisson distribution with the appropriate expectation value. Simple phase fluctuations can be applied to the fringes. Figure 31 shows two sweeps of example output.

c Main program for non-dispersed fake data
c Data-faking routine for COAST interferometer
c Written by Bob Tubbs, St John's College
c for Part III Physics Project, March 1998

include 'rnt.inc'

c Parameters

real pi
parameter(pi = 3.14159265358979)

c Variables

integer sample, subsample, sweep, nsweeps
integer count1, count2, direction

real spacing, position, filterwidth, intensity
real mean, amplitude, seed, initialfringeposition
real fringeposition, fringevelocity, fringeacceleration
real loggaussian, gaussian, time

character outfile*80, photonnoiseflag*8

c Function names

integer photons
real rand

c Start

status = ok

write(6,*) "Fringe spacing (samples)?"
read(5,*) spacing
write(6,*) "Initial position of fringe envelope centre",
* " (samples)?"
read(5,*) initialfringeposition
write(6,*) "Initial velocity of fringe envelope (samples/s)?"
read(5,*) fringevelocity
write(6,*) "Acceleration of fringe envelope (samples/s**2)?"
read(5,*) fringeacceleration
write(6,*) "Fractional bandwidth of filter?"
write(6,*) "(Between 1/e points)"
read(5,*) filterwidth
write(6,*) "Mean Intensity"
read(5,*) mean
write(6,*) "Maximum intensity fluctuation"
read(5,*) amplitude
write(6,*) "Number of sweeps"
read(5,*) nsweeps
write(6,*) "Output filename?"
read(5,*) outfile
write(6,*) "Add photon noise (y/n)?"
read(5,*) photonnoiseflag
if ((photonnoiseflag .eq. "y") .or.
* (photonnoiseflag .eq. "Y")) then
photonnoiseflag = "y"

c Seed psuedo-random number generator

write(6,*) "Random number seed?"
read(5,*) seed
seed = rand(seed)
endif

c Initial trolley direction

direction = 1

c Open output file

open(UNIT=9, STATUS='NEW', ERR=90, FILE=outfile)

c Output number of sweeps to file

write(9,FMT=*) nsweeps

c Start main loop

do sweep = 1, nsweeps

do sample = 1, 500

c Calculate mean light intensity for this sample

intensity = 0.0

do subsample = 1, 5

c Calculate time since first sample (seconds)
c There are 10 sweeps per second, 5000 samples
c per second, 25000 subsamples per second

time = real(sweep - 1)/10.0 + real(sample - 1)/5000.0
* + real(subsample - 3)/25000.0

c Calculate current fringe position (samples)

fringeposition = initialfringeposition + fringevelocity *
* time + fringeacceleration * (time**2) / 2.0

c Calculate current trolley position (samples)

position = 250.0 - 250.0 * real(direction) +
* (real(sample - 1) + real(subsample - 3)/5.0)
* * real(direction)

c Formula for log(Gaussian) centred at (position - fringeposition)

loggaussian = -((position - fringeposition) *
* filterwidth/spacing)**2

c Calculate Gaussian, checking for floating point underflow

if (loggaussian .gt. -30.0) then
gaussian = exp(loggaussian)
else
gaussian = 0.0
endif

c Calculate change in intensity produced by interference
c Add 1/5 of this value to the total value for the sample
c Repeated for each of the 5 subsamples to give mean value

intensity = intensity + amplitude*gaussian
* *cos(2*pi*(position - fringeposition)/spacing)/5.0

enddo

c Calculate photon counts for case with beams in phase and
c case for beams in anti-phase

if (photonnoiseflag .eq. "y") then

c Use Poisson distributed random number generator: photons()

count1 = photons(mean + intensity)
count2 = photons(mean - intensity)

else

count1 = int(mean + intensity)
count2 = int(mean - intensity)

endif

c Output photon counts to file

write(9,FMT=*) sample + (sweep - 1) * 500, count1, count2

enddo

c Reverse trolley direction

direction = - direction

c End main loop

enddo

c Close output file

close(UNIT=9, ERR=90)

stop

90 write(6,*) "Filer error"

end

 
 
c Subroutine which generates Poisson distributed
c pseudo-random numbers. Ignoring deficiencies in
c the FORTRAN random number generator and rounding
c errors, the psuedo-random numbers are perfectly
c Poisson-distributed.

integer function photons(intensity)

implicit none

c Variables

real intensity

integer i

real r, poisson, cumulative, logfactoriali, logpoisson

c Functions

real rand

c Generate psuedo-random number with rectangular distribution

r = rand(0.0)

c Initialise variables

i = 0
logfactoriali = 0.0
cumulative = 0.0
poisson = 1.0

c Check for intensity less than or equal to zero

if (intensity .gt. 0.0) then

c Main loop
c Repeated until cumulative probability is greater than r or has
c reached maximum value. Total cumulative probability may be
c minutely smaller than 1.0 due to rounding errors etc.

do while ((cumulative .lt. r) .and.
* ((real(i) .le. intensity) .or. (poisson .gt. 0.0)))

logpoisson = -intensity+real(i)*log(intensity)-logfactoriali

c Calculate probability of i photons being emmitted, checking for
c floating point underflow

if (logpoisson .gt. -30.0) then

poisson = exp(logpoisson)

else

poisson = 0.0

endif

c Calculate cumulative probability of between 0 and i photons being
c emmitted.

cumulative = cumulative + poisson

c Increment i and calculate new value of log(factorial(i))

i = i + 1

logfactoriali = logfactoriali + log(real(i))

enddo

c Number of photons required is i - 1

photons = i - 1

else

photons = 0

endif

return

end

Figure 31 - Example output from the Fake data routine. Two sweeps each of length 500 samples are shown. The mean intensity was 100 photons per sample. The peak fringe visibility was 0.5. The fractional bandwidth of the optical filter was 0.05.

Return to Bob's Interferometry Page.

Return to Bob's Home Page.