I am doing an atmospheric chemistry PhD and I use Fortran to model the reactions that I do in a flow tube.
I'm not really the best programmer, I can alter the model to my needs, but I wouldn't be able to write anything from scratch.
I use Monte Carlo simulations on my reactions in order to work out the errors in the rate coefficients. At the moment i use the code below and it does a "top hat" shaped distribution. Obviously there is more to model, as it actually models the fit to my data as well etc.
! here are parameters
If (zed.lt.0) then
write(9,*) 'N = ',N
write(*,*) 'error n = ',mcerrorn
write(9,*) 'error n = ',mcerrorn
write(9,*) 'sum n = ',Xn
write(*,*) 'error p = ',mcerrorp
write(9,*) 'error p = ',mcerrorp
write(9,*) 'sum p = ',Xp
Basically I work with rate coefficients and concentrations so there is a variable for each parameter in my model that has an error associated with it. The model calculates the overall error
However I need a normal distribution. I have been given some code by my supervisor (obviously not in Fortran...) but I have no idea about how to put it into my model.
17 x = mean + 25*sd*(0.5-RND)*2
p = exp(-(x-mean)^2/2/sd^2)
if RND > p the goto 17
What my supervisor said about the above code: I'm to assume the parameter x has a best fit value "mean" and is normally distributed with a standard deviation (or standard error ("sd").
In line 17 x is chosen to be somewhere within +/- 25 standard deviations of the mean (RND is a random number between 0 and 1). The probability p of x with respect to the mean is then calculated according to the normal distribution. A second random number is the generated. If RND < p, then the resulting value of x is used. Otherwise the routine loops back to line 17 and repeats.
Sorry if this makes sense of I haven't included enough information. Please feel free to ask me any questions.
Thanks in advance for any help