- #1
karenmarie3
- 10
- 0
Hey everyone,
I've been working on this project for awhile (some of you may recall having seen a very early version), and I am so close to being finished, but I have a problem somewhere, and I can't figure out why.
So what is happening here is that this code is supposed to calculate the average displacement of an object from some point chosen randomly. In the first steps of writing this code, I started with simply the choice of random points, and then a third random point was chosen and compared to the Gaussian calculated by the original 2 points. (First "do" loop.) When I did this part, everything worked fine. My points all had a higher probability of being created in the center of the circle I have created in the program.
The problem is, that now when I run this with all my extra steps, I am getting a large number of lengths that are short, instead of long... suggesting most of my creation points are no longer near the center, which is not the case.
I have been looking at this for way too long, and I cannot find the problem. This is the first code I have written completely on my own, in any language, so I may be missing subtleties of the language or something related to syntax. I realize there are some possibly more efficient ways of doing things than the way I am, but could someone please give me some direction as to what is going wrong, or where my problem is? I don't necessarily want to be told what to change as to where to look right now.
Thanks!
I've been working on this project for awhile (some of you may recall having seen a very early version), and I am so close to being finished, but I have a problem somewhere, and I can't figure out why.
So what is happening here is that this code is supposed to calculate the average displacement of an object from some point chosen randomly. In the first steps of writing this code, I started with simply the choice of random points, and then a third random point was chosen and compared to the Gaussian calculated by the original 2 points. (First "do" loop.) When I did this part, everything worked fine. My points all had a higher probability of being created in the center of the circle I have created in the program.
The problem is, that now when I run this with all my extra steps, I am getting a large number of lengths that are short, instead of long... suggesting most of my creation points are no longer near the center, which is not the case.
I have been looking at this for way too long, and I cannot find the problem. This is the first code I have written completely on my own, in any language, so I may be missing subtleties of the language or something related to syntax. I realize there are some possibly more efficient ways of doing things than the way I am, but could someone please give me some direction as to what is going wrong, or where my problem is? I don't necessarily want to be told what to change as to where to look right now.
Thanks!
Code:
program avglength
implicit none
integer N, L, M, i, j, k
double precision gauss_dist, a, a2, z_min, e, total
double precision z_max, time, tmprt, zlength, y_length
double precision avg_length, zcount, rad, x_length
double precision x, y, z, x_path, y_path, beta, pi
double precision delta_t, zlength_sum
open(unit=2, file="x-y-acc.dat", status="unknown")
open(unit=3, file="lengths.dat", status="unknown")
open(unit=4, file="x-y-z.dat", status="unknown")
open(unit=6, file="x-and-y_path_beta.dat", status="unknown")
open(unit=7, file="avg_length.dat", status="unknown")
open(unit=8, file="x_and_y_path.dat", status="unknown")
open(unit=9, file="x_length-y_length.dat", status="unknown")
open(unit=11, file="x-y-rad.dat", status="unknown")
open(unit=12, file="count.dat", status="unknown")
! value assignment for counters & constants !
N = 10000
M = 100
L = 135
a = 7
a2 = a*a
e = 2.718
pi = 3.1415
time = 0.0d0
delta_t = 0.2d0
zcount = 0.0d0
!
! !
! creation and motion evolution !
! !
!
! selection of seed points (seed loop) !
do i = 1, N
x = (rand() - 0.5) * 26.0
y = (rand() - 0.5) * 26.0
z = rand()
gauss_dist = e**(-(x*x + y*y)/a2)
write(4,*) x, y, z, gauss_dist
if (z .le. gauss_dist) then
write(2,*) x, y
end if
if (z .le. gauss_dist) then
! assignment of M angles to each jet seed (angle loop) !
do j = 1, M
time = 0.0d0
beta = 2.0*pi*rand() !angle of expulsion
x_path = x
y_path = y
write(6,*) x_path, y_path, beta
! time evolution of movement (time loop) !
do k = 1, L
rad = dsqrt(x_path*x_path + y_path*y_path)
write(11,*) x_path, y_path, rad
if ( rad .le. 10.0d0 ) then
tmprt = 1.0d0
else if (rad .gt. 10.0d0 ) then
tmprt = 0.0d0
end if
if ( tmprt .gt. 0.16) then
time = time + delta_t
x_path = x_path + delta_t*cos(beta)
y_path = y_path + delta_t*sin(beta)
write(8,*) x_path, y_path, time
else if (tmprt .le. 0.16) then
exit
! preceding two lines should cause computer to skip out of this loop to calculate jet length once temp drops below 0.16GeV
end if
end do
x_length = x_path - x
y_length = y_path - y
write(9,*) x_length, y_length
zlength = dsqrt (x_length**2 + y_length**2)
zcount = zcount + 1.00
write(3,*) zlength, zcount
end do
end if
end do
close(3)
close(4)
close(6)
close(8)
close(11)
!
! !
! summation of lengths followed by averaging (counting loop) !
! !
!
! initilaization of values in counting loop !
zlength_sum = 0.0d0
write(12,*) zcount
open(unit=3, file="lengths.dat", status="unknown")
N = 1000
do i = 1, N
read(3,*) zlength
zlength_sum = zlength_sum + zlength
end do
avg_length = zlength_sum/zcount
write(7,*) avg_length
stop
close(7)
close(12)
end program