 /*
  * Khoros: $Id: noise.c,v 1.1 1991/05/10 15:41:26 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: noise.c,v 1.1 1991/05/10 15:41:26 khoros Exp $";
#endif

 /*
  * $Log: noise.c,v $
 * Revision 1.1  1991/05/10  15:41:26  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, University of New Mexico.  All rights reserved.
 * 
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *----------------------------------------------------------------------
 */

#include "unmcopyright.h"	 /* Copyright 1990 by UNM */
#include "vinclude.h"	

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>	    file name: noise.c                                <<<<
   >>>>                                                       <<<<
   >>>>   description: contains routines that generate 2d     <<<<
   >>>>                noise "images"  If either the row or   <<<<
   >>>>                column arguments are set to 1, then    <<<<
   >>>>                a 1d noise data set is generated.      <<<<
   >>>>                In all cases, the data is stored in    <<<<
   >>>>                float vector that must be recast as    <<<<
   >>>>                necessary.                             <<<<
   >>>>                                                       <<<<
   >>>>                This is part of vmath.a                <<<<
   >>>>                                                       <<<<
   >>>>      routines: ggauss(), gexpon(), gunif(), grleigh() <<<<
   >>>>                gpoiss()  			      <<<<
   >>>>                                                       <<<<
   >>>> modifications: See module headers         	      <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

#define LN10 2.3025851


/***********************************************************************
*
*  Routine Name: ggauss()
*
*          Date: April 23, 1990
*        
*       Purpose: Generates a 2d Gaussian random noise data set 
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                mean     - (float) mean of the data set
*                variance - (float) variance of the data set 
*
*        Output: 1 - error condition
*                0 - no error 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int ggauss(vect,rows,cols,mean,variance)
    float *vect,mean,variance;
    int rows,cols;

{
    int length,i;
    float temp,u1,u2,uno,urng();
    
    length = rows * cols;

    uno = 2.0 * exp(0.825);

   /* generate the noise sequence */
        for(i = 0; i < length/2; i++)
        {
    	  u1 = urng(); 
  	  u2 = urng();
          temp = sqrt( (double) (-uno*variance*log10( (double)u1 )));
	  vect[2*i] = temp * cos(XV_2PI*u2) + mean;
	  vect[2*i+1] = temp * sin(XV_2PI*u2) + mean;
        }
 
    if (length & 1) {
           u1 = urng();
           u2 = urng();
           temp = sqrt( (double) (-uno*variance*log10( (double)u1 )));
           vect[length-1] = temp * cos(XV_2PI*u2) + mean;
    }

    return(1);
}

/***********************************************************************
*
*  Routine Name: gexpon()
*
*          Date: May 10, 1990
*        
*       Purpose: Generates a 2d exponential random noise data set 
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                variance - (float) variance of the data set 
*
*        Output: 1 - error condition
*                0 - no error 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int gexpon(vect,rows,cols,variance)
    float *vect,variance;
    int rows,cols;

{
    int length,i;
    float u1,urng();
    
    length = rows * cols;

/*
** generate the noise sequence
*/

        for(i = 0; i < length; i++)
        {
    	  u1 = urng(); 
	  vect[i] = -1.0*sqrt(variance)*LN10*log10(u1);
        }
 
    return(1);
}

/***********************************************************************
*
*  Routine Name: gunif()
*
*          Date: April 24, 1990
*        
*       Purpose: Generates a 2d uniform random noise data set 
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                minimum  - (float) minimum of the data set
*                maximum  - (float) maximum of the data set 
*
*        Output: 1 - error condition
*                0 - no error 
*
*    Written By: Jeremy Worley 
*
* Modifications: Fri Aug 31 15:04:34 MDT 1990 Jeremy Worley
*		 An awful bug that i'm quite embarrassed about:  original
*		 amplitude calculation read "amp = minimum-maximum".
*
***********************************************************************/

int gunif(vect,rows,cols,minimum,maximum)
    float *vect,minimum,maximum;
    int rows,cols;

{
    int length,i;
    float u1,amp,urng();

    length = rows * cols;
    amp = maximum - minimum;

        for(i = 0; i < length; i++)
        {
          u1 = urng();
          vect[i] = amp * u1 + minimum;
        }

    return(1);
}

/***********************************************************************
*
*  Routine Name: grleigh()
*
*          Date: May 2, 1990
*        
*       Purpose: Generates a 2d Rayleigh random noise data set 
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                variance - (float) variance of the data set 
*
*        Output: 1 - error condition
*                0 - no error 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int grleigh(vect,rows,cols,variance)
    float *vect,variance;
    int rows,cols;

{
    int junk1,junk2,length,i,ggauss();
    float *t1,*t2;
        
    length = rows * cols;

/*
** allocate some space for gaussian temporary random variables
*/

    t1 = (float *)malloc(length*sizeof(float));
    if(t1==NULL){
       fprintf(stderr,"grleigh(math.a):  unable to malloc temp space 1.\n");
       return(0);
    }
    t2 = (float *)malloc(length*sizeof(float));
    if(t2==NULL){
       fprintf(stderr,"grleigh(math.a):  unable to malloc temp space 2.\n");
       free(t1);
       return(0);
    }

/*
** generate the noise sequence
*/

    junk1 = ggauss(t1,rows,cols,(float)0,(float)(variance/(2.0 - XV_2PI/4)));
    junk2 = ggauss(t2,rows,cols,(float)0,(float)(variance/(2.0 - XV_2PI/4)));
    if(!junk1 || !junk2){
       free(t1);
       free(t2);
       fprintf(stderr,"grleigh(math.a):  an error occurred during ggauss.\n");
       return(0);
    }

/*
** rayleigh is the square root of the sum of the squares of two gaussian
** random variables
*/

    for(i = 0; i < length; i++) {
          vect[i] = sqrt((double)(t1[i]*t1[i] + t2[i]*t2[i]));
    }

/*
** free up temporary space
*/

    free(t1);
    free(t2);
 
    return(1);
}


/***********************************************************************
*
*  Routine Name: gpoiss()
*
*          Date: May 21, 1990
*        
*       Purpose: Generates a 2d Poisson random noise data set 
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                variance - (float) variance of the data set 
*                time     - (float) time argument 
*
*        Output: 1 - error condition
*                0 - no error 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int gpoiss(vect,rows,cols,variance,time)
    float *vect,time,variance;
    int rows,cols;

{
    int length,i;
    float t,tsum,p,u1,urng();

    length = rows * cols;

/*
** generate the poisson random variable
*/

    for(i=0;i<length;i++){
        tsum = p = 0;
        do{
           u1 = urng();
           t  = (-1.0*time/variance)*LN10*log10(1.0 - u1);
           tsum += t;
           if(tsum <= time) p += 1.0;
        }while(tsum <= time);
       *vect++ = p;
    }

/*
** normal return
*/
    
    return(1);
}


