Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran] Filling a disk with random points

  1. Jun 12, 2014 #1
    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. jcsd
  3. Jun 12, 2014 #2

    jedishrfu

    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.
     
  4. Jun 12, 2014 #3

    mathman

    User Avatar
    Science Advisor
    Gold Member

    If you want the points uniform within the ring, choose r2 uniformly and take square root, rather than choosing r uniformly.
     
  5. Jun 12, 2014 #4
    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
     
  6. Jun 12, 2014 #5

    jedishrfu

    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.
     
  7. Jun 13, 2014 #6

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  8. Jun 14, 2014 #7

    mathman

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  9. Jun 14, 2014 #8
    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
     
  10. Jun 14, 2014 #9

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  11. Jun 15, 2014 #10
    Thanks FactChecker.
    I don't understand why my call to seed produce 0 every time I run the program?
     
    Last edited: Jun 15, 2014
  12. Jun 15, 2014 #11

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  13. Jun 15, 2014 #12
    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
     
  14. Jun 15, 2014 #13

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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
  15. Jun 15, 2014 #14
    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.
     
  16. Jun 15, 2014 #15

    jtbell

    User Avatar

    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.
     
  17. Jun 15, 2014 #16
    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?
     
  18. Jun 15, 2014 #17

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  19. Jun 15, 2014 #18

    mathman

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  20. Jun 15, 2014 #19

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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
     
     
  21. Jun 16, 2014 #20
    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




Similar Discussions: [Fortran] Filling a disk with random points
Loading...