:::::::::::::: README :::::::::::::: This program will produce a raster map layer of gaussian deviates whose mean and standard deviation can be expressed by the user. It is essentialy the same as r.rand.surf, but uses a gaussian random number generator instead. Both random number generators are taken from Press, Flannery, Teukolsky and Vetterling (1988) - Numerical Recipes in C. [JWO 24.10.91] :::::::::::::: main.c :::::::::::::: /* ** Code Compiled by Jo Wood [JWO] 24th October 1991 ** Midlands Regional Research Laboratory (ASSIST) ** ** */ #include "gis.h" main(argc,argv) int argc; char *argv[]; { /****** INITIALISE ******/ double gauss_mean,gauss_sigma; struct Option *out; /* Structures required for the G_parser() */ /* call. These can be filled with the */ struct Option *mean; /* various defaults, mandatory paramters */ /* etc. for the GRASS user interface. */ struct Option *sigma; G_gisinit (argv[0]); /* This GRASS library function MUST be called first to check for valid database and mapset. As usual argv[0] is the program name. This can be recalled using G_program_name(). */ /****** SET PARSER OPTIONS ******/ out = G_define_option(); /* Request pointer to memory for each option. */ mean = G_define_option(); /* Mean of the distribution */ sigma = G_define_option(); /* Standard deviation of the distribution */ out->key = "out"; out->description = "Name of the random surface to be produced"; out->type = TYPE_STRING; out->required = YES; mean->key = "mean"; mean->description = "Distribution mean"; mean->type = TYPE_DOUBLE; mean->answer = "0.0"; sigma->key = "sigma"; sigma->description = "Standard deviation"; sigma->type = TYPE_DOUBLE; sigma->answer = "1.0"; if (G_parser(argc,argv)) exit(-1); /* Returns a 0 if sucessful */ sscanf(mean->answer,"%lf",&gauss_mean); sscanf(sigma->answer,"%lf",&gauss_sigma); /****** CHECK THE CELL FILE (OUT) DOES NOT ALREADY EXIST******/ if (G_legal_filename(out->answer)==NULL) { char err[256]; sprintf(err,"Illegal file name. Please try another."); G_fatal_error(err); } else { if (G_find_cell(out->answer,"") !=NULL) { char err[256]; sprintf(err,"Raster map [%s] already exists.\nPlease try another.",out); G_fatal_error(err); } } /****** CREATE THE RANDOM CELL FILE ******/ gaussurf(out->answer,gauss_mean,gauss_sigma); } :::::::::::::: gaussurf.c :::::::::::::: /*****************/ /** gaussurf() **/ /*****************/ #include "gis.h" #include gaussurf(out,mean,sigma) char *out; /* Name of cell files to be opened. */ double mean,sigma; /* Distribution parameters. */ { int nrows,ncols; /* Number of cell rows and columns */ CELL *row_out; /* Buffer just large enough to hold one */ /* row of the raster map layer. */ int fd_out; /* File descriptor - used to identify */ /* open cell files. */ int row_count,col_count; float rand1(),gauss(); /****** INITIALISE RANDOM NUMBER GENERATOR ******/ rand1(-1* getpid()); /****** OPEN CELL FILES AND GET CELL DETAILS ******/ fd_out = G_open_cell_new(out); nrows = G_window_rows(); ncols = G_window_cols(); row_out = G_allocate_cell_buf(); /****** PASS THROUGH EACH CELL ASSIGNING RANDOM VALUE ******/ for (row_count=0;row_count float gauss(seed) int seed; { static int iset=0; static float gset; float fac,r,v1,v2; float rand1(); if (iset==0) { do { v1=2.0*rand1(seed)-1.0; v2=2.0*rand1(seed)-1.0; r=v1*v1+v2*v2; } while (r>=1.0); fac=sqrt(-2.0*log(r)/r); gset=v1*fac; iset=1; return(v2*fac); } else { iset=0; return gset; } } :::::::::::::: rand1.c :::::::::::::: /****************************************************************/ /*** Random Number Generator (Uniform Deviates 0.0 -> 1.0) ***/ /*** ***/ /*** Based on three linear congruential generators: ***/ /*** One for most significant part, ***/ /*** One for the least significant, ***/ /*** One for a shuffling routine (Knuth, 1981). ***/ /*** [Press et al, 1988] ***/ /*** ***/ /*** Coded Oct 23 1991 ***/ /*** Version 1.0 ***/ /*** ***/ /****************************************************************/ /* Arbitrary constants */ #define M1 259200 #define IA1 7141 #define IC1 54773 #define RM1 (1.0/M1) #define M2 134456 #define IA2 8121 #define IC2 28411 #define RM2 (1.0/M2) #define M3 243000 #define IA3 4561 #define IC3 51349 float rand1(seed) int seed; { static long ix1,ix2,ix3; float temp; int j; static float r[98]; static int iff=0; if (seed < 0 || iff==0) /* Initialise if negative seed is given. */ /* This MUST be done on the first call. */ { iff=1; ix1=(IC1-seed) % M1; /* Seed the 1st routine */ ix1=(IA1*ix1+IC1) % M1; /* and use it to seed the */ ix2=ix1 % M2; /* 2nd and 3rd routines */ ix1=(IA1*ix1+IC1) % M1; ix3=ix1 % M3; for (j=1; j<=97; j++) { ix1=(IA1*ix1+IC1) % M1; /* Fill table with sequen- */ ix2=(IA2*ix2+IC2) % M2; /* tial uniform deviates */ r[j]=(ix1+ix2*RM2)*RM1; /* from the 1st 2 routines */ } seed=1; } ix1=(IA1*ix1+IC1) % M1; /* Generate the next number in the */ ix2=(IA2*ix2+IC2) % M2; /* sequence of routine */ ix3=(IA3*ix3+IC3) % M3; j=1 + ((97*ix3) / M3); /* Get random position ion array */ temp = r[j]; r[j]=(ix1+ix2*RM2)*RM1; /* Refil array and return entry */ return(temp); }