[Fortran] Filling a disk with random points

[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

 

[jedishrfu]

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.

 

[mathman]

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

 

[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

 

[jedishrfu]

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.

 

[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.

 

[mathman]

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

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.

 

[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

 

[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.

 

[mrcotton]

Thanks FactChecker.
I don't understand why my call to seed produce 0 every time I run the program?

 

[FactChecker]

Thanks FactChecker.
I don't understand why my call to seed produce 0 every time I run the program?

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.

 

[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

 

[FactChecker]

Hi FactChecker,
So the line in my code
call srand(seed) does not generate the value of seed then?
what does it do?
D

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.

 

[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.

 

[jtbell]

In general if I see subroutine(value) then value is passed to the sub routine?

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.

 

[mrcotton]

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

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?

 

[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.

 

[mathman]

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?

I am not sure what you said. My interpretation:
Let [a,b] be the range for r. Calculate x uniformly in [a²,b²], then r = √x.

 

[FactChecker]

I wrote a short program to test this on 100 million points with radius between 0.5 and 2.0

$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.

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

 

[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?

 

[mrcotton]

Hi Factchecker, Thanks for your time on this. I can see how this is a good way to check if my generated radii are evenly distributed. Please excuse my ignorance but what language is this?
Is there unseen code that generates the count for each bin (if that is the correct word) for the .1 increments. What exactly is the scaled count? and the scaled 1/r?
Thanks again
D

 

[mrcotton]

The story so far,
I must thank you all I could not have got this far without you guys!
I need to work out how to plot my points in gnu plot next and how to store my points in an array

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

real,parameter :: pi = 3.1415926
real,parameter :: r_max = 10
real,parameter :: r_min = 1

integer :: seed, i
! radius rad are polar coordinates angle in radians
! x,y are cartesian coordinates

real :: radius, rad, x, y
seed = 7873
i=1

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

DO i = 1, 100
! calculate the radius
radius = (r_max-r_min)*rand()
! calculate the polar angle in radians
rad = 2*pi*rand()

x=radius*cos(rad)
y=radius*sin(rad)

open (unit = 2, file = "coordinates.txt")
Write (2,*) x,y

write (*,*) 'radius',radius
write (*,*) 'polar',rad
write (*,*) 'x',x
write (*,*) 'y',y
END DO


end program test_rand

 

[mathman]

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?

The area differential in polar coordinates is dA = rdrdθ.

 

[jtbell]

Please excuse my ignorance but what language is this?

It's Perl.

Is there unseen code that generates the count for each bin (if that is the correct word) for the .1 increments.

It looks like he put some of the code for the histogram into a subroutine named &record_hist, and omitted a bunch of stuff after the code that he does show.

 

[FactChecker]

Hi Factchecker, Thanks for your time on this. I can see how this is a good way to check if my generated radii are evenly distributed. Please excuse my ignorance but what language is this?
Is there unseen code that generates the count for each bin (if that is the correct word) for the .1 increments. What exactly is the scaled count? and the scaled 1/r?
Thanks again
D

@jtbell is right. I didn't want the bookkeeping to distract from the essential point. Here is the omitted code at the end of the program:

sub record_hist {
	# record results in a histogram for every 0.1 radius increment
	$hist_index = int(10*$rand_radius);
	$hist[$hist_index]++;
}

# print the histogram results between the original radius limits
print <<END;
Histogram results
  min r    max r    count      Scaled  Scaled ~ 1/r
                                        count
END
foreach $i (5..19){
	# histogram cell range of radius
	$hist_lower = $i/10;
	$hist_upper = ($i+1)/10;

	# scale raw numbers versus cell for least radius (should grow proportional to radius)
	$scaled_to_first_cell = $hist[$i]/$hist[5];

	# scale again versus inverse of middle radius of the cell (should be 1.0)
	$scaled_to_inverse_radius = $scaled_to_first_cell * ( 0.55/($hist_lower+0.05));

	printf "%6.3f  %6.3f  %8d  %6.3f  %6.3f\n",
	   $hist_lower,
	   $hist_upper,
	   $hist[$i],
	   $scaled_to_first_cell,
	   $scaled_to_inverse_radius;
}

 

[jedishrfu]

Not to be left out, here's a Processing IDE sketch using polar coordinates to randomly populate the annulus:

// DEFINE canvas size
int xsize=400;
int ysize=400;

// COMPUTE center of canvas
int xhalf=xsize/2;
int yhalf=ysize/2;

// DEFINE inner and outer radii
float rinner=50;
float router=200;

// DEFINE number of random dots to draw
int NDOTS=100;
int dots=0;

void setup() {
  // DEFINE a xsize pixels x ysize pixels canvas
  size(xsize,ysize);
  
  // DRAW outer circle with diameter 2*router pixels
  ellipse (xhalf,yhalf,2*router,2*router);
  
  // DRAW inner circle with diameter 2*rinner pixels and fill in as light GRAY
  fill(200,200,200);
  ellipse(xhalf,yhalf,2*rinner,2*rinner);
}

void draw() {
 
  // DRAW no more than NDOTS
  if(dots<NDOTS) dots++; else return;
  
  // COMPUTE polar position of random dot
  float r = (float)(randomWithinRange(rinner,router));
  float theta=(float)(randomWithinRange(0.0,2*Math.PI));
  
  //println("r="+r+"  theta="+theta);
  
  // COMPUTE x,y position of random dot
  float xdot=r*cos(theta)+xhalf;
  float ydot=r*sin(theta)+yhalf;
  
  // DRAW a random dot 2 pixel diameter
  ellipse(xdot,ydot,2,2);

}

// ===============================================================
// The following method: randomWithinRange(min,max) was copied from this webpage:
//
//         http://stackoverflow.com/questions/7961788/math-random-explained
// ===============================================================

// RETURN a random value between two values
double randomWithinRange(double min, double max)
{
   double range = Math.abs(max - min);     
   return (Math.random() * range) + (min <= max ? min : max);
}

To run the sketch download the IDE from Processing.org, its designed for casual programmers and graphics artists of all ages. :-)

 

[mrcotton]

Thanks guys, of all ages eh, including senile old codgers like me. There is hope for this old dog yet!