[Dart-dev] [6963] DART/trunk/random_seq: change the format of the output in test_diff to be

[nancy at ucar.edu]

Revision: 6963
Author:   nancy
Date:     2014-04-30 15:29:06 -0600 (Wed, 30 Apr 2014)
Log Message:
-----------
change the format of the output in test_diff to be
more clear about what's computed and what the inputs are.
change test_hist so it writes the binned output into
a matlab script instead of printing to the screen.

Modified Paths:
--------------
    DART/trunk/random_seq/test_diff.f90
    DART/trunk/random_seq/test_hist.f90

-------------- next part --------------
Modified: DART/trunk/random_seq/test_diff.f90
===================================================================
--- DART/trunk/random_seq/test_diff.f90	2014-04-30 21:24:12 UTC (rev 6962)
+++ DART/trunk/random_seq/test_diff.f90	2014-04-30 21:29:06 UTC (rev 6963)
@@ -30,6 +30,9 @@
 type (random_seq_type) :: r
 
 real(r8) :: r1, r2, dist2, mean_dist2
+character(len=16) :: wformatI = '(A,I8)'
+character(len=16) :: wformat1 = '(A,1(F25.16))'
+character(len=16) :: wformat2 = '(A,2(F25.16))'
 
 !-----------------------------------------------------------------------------
 ! Namelist with default values
@@ -47,23 +50,25 @@
 read(iunit, nml = test_diff_nml, iostat = io)
 call check_namelist_read(iunit, io, "test_diff_nml")
 
-write(*, *) 'sample size = ', n
-write(*, *) 'mean 1, 2 = ', 0.0, mean
-write(*, *) 'stddev 1, 2  = ', sd1, sd2
-write(*, *) 'predicted distance = ', sqrt(mean**2 + sd1**2 + sd2**2)
 write(*, *) ''
+write(*, wformatI) 'sample size = ', n
+write(*, *) ''
+write(*, wformat2) 'mean, stddev 1 = ', 0.0_r8, sd1
+write(*, wformat2) 'mean, stddev 2 = ', mean, sd2
+write(*, *) ''
 
 
 call init_random_seq(r,iseed)
-mean_dist2 = 0.0d0
+mean_dist2 = 0.0_r8
 do i = 1, n
-   r1 = random_gaussian(r, 0.0d0, sd1)
-   r2 = random_gaussian(r,  mean, sd2)
+   r1 = random_gaussian(r, 0.0_r8, sd1)
+   r2 = random_gaussian(r,  mean,  sd2)
    dist2 = (r1 - r2)**2
    mean_dist2 = mean_dist2 + dist2
 end do
 
-write(*, *) 'sample mean distance = ', sqrt(mean_dist2 / n)
+write(*, wformat1) 'predicted distance   = ', sqrt(mean**2 + sd1**2 + sd2**2)
+write(*, wformat1) 'sample mean distance = ', sqrt(mean_dist2 / n)
 write(*, *) ''
 
 call finalize_utilities()

Modified: DART/trunk/random_seq/test_hist.f90
===================================================================
--- DART/trunk/random_seq/test_hist.f90	2014-04-30 21:24:12 UTC (rev 6962)
+++ DART/trunk/random_seq/test_hist.f90	2014-04-30 21:29:06 UTC (rev 6963)
@@ -11,7 +11,8 @@
 
 use      types_mod, only : r4, r8, digits12
 use  utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
-                           initialize_utilities, finalize_utilities
+                           initialize_utilities, finalize_utilities,     &
+                           open_file, close_file
 use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform
 
 implicit none
@@ -24,7 +25,7 @@
 
 type (random_seq_type) :: r
 real(r8) :: val
-integer :: i, bin
+integer :: i, bin, iunit
 integer, parameter :: nbins = 500
 integer :: repcount = 100000000
 integer :: bincount(nbins)
@@ -32,22 +33,41 @@
 call initialize_utilities('test_hist')
 call register_module(source,revision,revdate)
 
+write(*, *) 'creating file "makehist.m" with output in matlab format'
+
+iunit = open_file("makehist.m", 'formatted', 'write')
+write(iunit, '(A)') 'bindata = [ ... '
 bincount(:) = 0
 
 call init_random_seq(r, 5)
 do i=1, repcount
    val = random_uniform(r)
-   bin = int(val*nbins)
+   ! generates a bin number between 0 and nbins-1, 
+   ! or possibly equal to nbins w/ roundoff error.  
+   ! (tests found 3 cases in 100M samples)
+   bin = floor(val*nbins) + 1
    !print *, i, val, bin
+   if (bin == nbins+1) bin = nbins
+   if (bin < 1 .or. bin > nbins+1) then
+      print *, 'error: computed bin = ', bin, ' should be > 1 and  < ', nbins
+   endif
    bincount(bin) = bincount(bin)+1
 enddo
 
 do i=1, nbins
    !print *, i, bincount(i)
    !print *, (1.0_r8/nbins)*(i-1), bincount(i)
-   print *, bincount(i)
+   write(iunit, '(I8,A)') bincount(i), ', ... '
 enddo
 
+write(iunit, '(A)') '];'
+write(iunit, '(A)') 'bar(bindata);'
+write(iunit, '(A)') 'axis([0, 500, 190000, 210000]);'
+
+write(*, *) 'closing "makehist.m" file'
+
+call close_file(iunit)
+
 call finalize_utilities()
 
 end program test_hist