# [Fortran] Filling a disk with random points

1. Jun 12, 2014

### mrcotton

Hi

My FORTRAN is rusty and my brain is even more rusty.

I want to populate an annulus with a randomly distributed set of points.

Any hints or tips to get me thinking about this would be gratefully received.

Thanks

D

2. Jun 12, 2014

### Staff: Mentor

Do you want them evenly spaced like a grid or more radially spaced?

Since its a disk then polar coordinates comes to mind. You could use the FORTRAN rand() function to get a random value from R-INNER to R-OUTER for the radial and use it again to get an angle from 0 to 2PI. Polar coordinates insures that your points fall within the boundary of the disk or ring figure.

3. Jun 12, 2014

### mathman

If you want the points uniform within the ring, choose r2 uniformly and take square root, rather than choosing r uniformly.

4. Jun 12, 2014

### mrcotton

Hi
Thanks to both of you for the advice.
Its going to be a steep learning curve.
I just re-partitioned my hard drive and have put on Ubuntu and am no looking for a FORTRAN 95 compiler.

I understand what you are saying about the polar coordinates and guess my next task is to work out how to use the rand() function to generate numbers within a range.

So it is better to calculate the square of r and then take the root. Excuse my ignorance but would this not just generate the same number?

Thanks for the help so far
D

5. Jun 12, 2014

### Staff: Mentor

for ubuntu:

http://stackoverflow.com/questions/5511704/how-can-i-install-fortran95-in-ubuntu-10-04

for the rand()

http://gcc.gnu.org/onlinedocs/gfortran/RAND.html

it returns a number from 0 to .99999999 but never 1

for a number from 0 to 99 then rand()*100 and assign it to an integer

notice there is the srand for seeding it with a value to generate a reproduceable random sequence for testing and then later you switch to full randomness.

http://infohost.nmt.edu/tcc/help/lang/fortran/rand.html

lastly they actually suggest that you use the RANDOM_NUMBER function as it produces a better more random number. Poor random number generators will have some sort of weakness like if you take every three numbers and plot them in 3-space then they will appear in the same plane...

Next do you have to use fortran if not then consider python or java. Python can handle large numbers and has been used in some scientific numeric applications. Python is readily available on Ubuntu and doesn't need to be installed. Processing.org is a good way to explore java interactively using graphics and comes with lots of examples.

And then there's matlab or one of its free clones like octave or freemat that could easily knock off this problem in a few lines of code.

http://freemat.sourceforge.net/

I used an earlier version of fortran called fortran-y which is roughly similar to fortran-77 I think.

6. Jun 13, 2014

### FactChecker

1) Decide if you want a) the polar coordinates or b) the cartesian coordinates to be uniformly spread.
a) Polar coordinates: angle = 2*pi*rand() and radius = (r_max-r_min)*rand().
b) Cartesian coordinates: x = (x_max-x_min)*rand() and y = (y_max-y_min)*rand(). Keep trying x and y till you get one that has an acceptable radius. --- That is a keeper. Ignore the others.

Use a uniform (0 to 1) random number generator, rand(), that satisfies your requirements. Some are better than others if your requirements are strict.

7. Jun 14, 2014

### mathman

The numbers will have the same range in either case. However the distribution will be different. If you generate the squares uniformly, the distribution of the points will be uniform in the ring. In other words you will get more points further from the center.

8. Jun 14, 2014

### mrcotton

I am very grateful to you all for the help so far.
It has been over 10 years since I wrote any code,
and as expected it runs unexpected.
Is this idea of the radius^2 that I generate r^2 then take the sqrt for my value of the radius?

+++++++++++++++++++++++++++++++++++++++++++++++++
program test_rand

real,parameter :: pi = 3.1415926
real,parameter :: r_max = 10
real,parameter :: r_min = 1
integer :: seed
real :: polar

call srand(seed)
write (*,*) 'seed',seed

! calculate the radius
radius = (r_max-r_min)*rand()
! calculate the radius
polar = 2*pi*rand()
write (*,*) 'radius',radius
write (*,*) 'polar',polar

end program test_rand
++++++++++++++++++++++++++++++++++++++++++++++++++

I have tried to read the links ablout rand and seeding.
I am as dumb as I look I am afraid.
I am missing the point somewhere as my code generates the same numbers each time.

**************************************************************
fortran@fortran-Latitude-E5500:~/FORTRAN/rand$./a.out seed 0 radius 2.1832044 polar 8.46310928E-02 fortran@fortran-Latitude-E5500:~/FORTRAN/rand$ ./a.out
seed 0
radius 2.1832044
polar 8.46310928E-02
fortran@fortran-Latitude-E5500:~/FORTRAN/rand$************************************************************** Do I need to put rand(seed) ? Thanks again D 9. Jun 14, 2014 ### FactChecker The seed is the starting point for the generation of a random series. You only set it once, before the loop starts then you let the rand function do whatever it wants with it inside the loop. Do not interfere with it after the first call. Given a seed number, the generated series will be identical every time you run the program. Not all seeds are good. Some functions need an initial integer seed that is not even. Google or the documentation might suggest a good seed. On another issue, you can generate uniform random points for any shape in Cartesian X-Y easily by generating a random X and a random Y so that the range of (X,Y) covers the area. Then reject any point that is not in the area of interest. The loss of efficiency due to rejecting random points is usually negligible. 10. Jun 15, 2014 ### mrcotton Thanks FactChecker. I don't understand why my call to seed produce 0 every time I run the program? Last edited: Jun 15, 2014 11. Jun 15, 2014 ### FactChecker Each seed number will produce the same series of random numbers on every run. That allows you to compare two things on two different runs using exactly the same random numbers. Seed it only once at the beginning of a run and then each future call (without a seed) should generate the next number in the series for that seed. Zero is probably a very bad seed. The numbers generated may not pass some tests of randomness. My guess is that a large prime number is a good seed but a lot depends on the exact algorithm used. You do not have to seed it and it will use the program start time to make a seed. But then each run will have different random numbers and not be repeatable. 12. Jun 15, 2014 ### mrcotton Hi FactChecker, So the line in my code call srand(seed) does not generate the value of seed then? what does it do? Do I just need to input the value seed into my program maybe with a parameter statement May I also ask about something else I read about srand() SEED Shall be a scalar INTEGER(kind=4). I have never heard of kind. This came from the GNU page https://gcc.gnu.org/onlinedocs/gfortran/SRAND.html#SRAND I appear to have become obsessed with my code instead of my physics Thanks D 13. Jun 15, 2014 ### FactChecker When you call srand(seed), you are giving it the seed. Because you have not set seed to another value, seed is 0 in every run, and your call is always the same as seed(0). So your following calls to rand gave you the random number series associated with a seed of 0. Try something like "call srand(7873)". 7873 is just a prime number I grabbed without knowing if it is a good one. It's probably ok. If you want a different series of random numbers for each run, I think you can skip the call to srand and it will use the program start time to generate it's own seed and the series generated will be different each time. Last edited: Jun 15, 2014 14. Jun 15, 2014 ### mrcotton Brilliant, for all this information on the net sometimes it takes a human mind to assist another human mind eh! I get what is going on now. In general if I see subroutine(value) then value is passed to the sub routine? Thanks for such rapid responses. 15. Jun 15, 2014 ### jtbell ### Staff: Mentor No, sometimes the value is returned from the subroutine to the program (or other subroutine) that calls it. And sometimes both things happen: a value is passed to the subroutine, which then changes it and returns the new value. To find out exactly what happens, you need to look at the documentation for the subroutine, or read the subroutine's code if you can. 16. Jun 15, 2014 ### mrcotton Does this mean I would calculate my radius using rand() as though it was the range I want squared then square root the value to make a more evenly distributed range? 17. Jun 15, 2014 ### FactChecker Yes. That will give random radii whose density is proportional to the radius. So the density of points will be independent of radius when they are spread in a 2*pi*r circle. PS. - I tried a little experiment to verify. 18. Jun 15, 2014 ### mathman I am not sure what you said. My interpretation: Let [a,b] be the range for r. Calculate x uniformly in [a2,b2], then r = √x. 19. Jun 15, 2014 ### FactChecker I wrote a short program to test this on 100 million points with radius between 0.5 and 2.0 Code (Text):$r_min = 0.5;
$r_max = 2;$r_min2=$r_min**2;$r_max2=$r_max**2; foreach$i (0..100000000){
# calculate uniform random between the squares of the upper and lower limits
$unif = ($r_max2 -$r_min2)*rand() +$r_min2;
# take the square root for the final random radious
$rand_radius =$unif**0.5;
# record results in a histogram for every 0.1 radius increment
&record_hist;
}
... record, scale and print results

Here are the results in a histogram with cells for each 0.1 radius delta. The density is proportional to radius (last column is all ~1.0), which is what was desired since the points will be spread out in a circle.
Code (Text):

Histogram results
min r   max r   count     Scaled  Scaled ~ 1/r
count
0.500   0.600   2934512   1.000   1.000
0.600   0.700   3464278   1.181   0.999
0.700   0.800   3999929   1.363   1.000
0.800   0.900   4532096   1.544   0.999
0.900   1.000   5068879   1.727   1.000
1.000   1.100   5601693   1.909   1.000
1.100   1.200   6133230   2.090   1.000
1.200   1.300   6667574   2.272   1.000
1.300   1.400   7204386   2.455   1.000
1.400   1.500   7727856   2.633   0.999
1.500   1.600   8267815   2.817   1.000
1.600   1.700   8804165   3.000   1.000
1.700   1.800   9333855   3.181   1.000
1.800   1.900   9864373   3.362   0.999
1.900   2.000  10395360   3.542   0.999

20. Jun 16, 2014

### mrcotton

Thanks mathman, that is what I was trying to say, It seems so much clearer in mathematical symbols. Is it possible to explain why this gives a better uniformity?

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook