Program hcorrnoise

!	reads the ggamma matrix computed by spacestats
!	and computes correlated noise.
 
		USE MSIMSLMS
		USE PORTLIB

		Implicit None
		

		Integer Ichan,Jchan,k,ios,instr
		Integer nchan /20/
		
		Integer Iseed /72338/ ! just a number I made up.
		Real cnoise(20,1000000)
		real ggamma(20,20)

		Real*8 s(20),ss(20),Numerator
		Real*8 sxy(20,20)

		Real Cov(20,20)
		Real r(20,20)

		Integer Npts
		Real*8 elapsed_time

		write(6,*) 'enter npts'
		read(5,*) npts

		do ichan = 1 , 20
	 	 s(ichan) = 0.0
		ss(ichan) = 0.0
		 do jchan = 1 , 20
		  sxy(ichan,jchan) = 0.0
	 	 enddo
		enddo

	    
		write(6,*) '11 or 12'
		read(5,*) instr
		if(instr .eq. 11) then
		 open(3,file='c:\osse\h11ggamma.dat',status='old',form='formatted')    
		else
		 open(3,file='c:\osse\h12ggamma.dat',status='old',form='formatted')
		endif

		do ichan = 1 , nchan		
		 read(3,6105,iostat=ios,err=310) ichan,(ggamma(ichan,jchan),jchan=1,10)		 
		 read(3,6105,iostat=ios,err=310) ichan,(ggamma(ichan,jchan),jchan=11,nchan)
		enddo
6105    format(i3,10e13.5)

	    Write(6,*) 'NOAA ',Instr, ' NPTS ',Npts
		Write(2,*) 'NOAA ',Instr, ' NPTS ',Npts
		
		elapsed_time = timef()
		
		Do k = 1 , npts

		 Call Correrror(Iseed,Nchan,Ggamma,CNoise(1,k))

		 do ichan=1,20
		   s(ichan) = s(ichan) + cnoise(ichan,k)
          ss(ichan) =ss(ichan) + cnoise(ichan,k)*cnoise(ichan,k)
		  do jchan = 1 , 20
		   sxy(ichan,jchan) = sxy(ichan,jchan) + cnoise(ichan,k)*cnoise(jchan,k)
		  enddo
		 enddo

        enddo ! k
		
		elapsed_time = timef()
		write(6,*) 'elapsed time ', elapsed_time
		pause


		do ichan = 1 , 20
		 do jchan = 1 , 20
		  Numerator=(float(npts)*sxy(ichan,jchan)-s(ichan)*s(jchan))
		  Cov(ichan,jchan) = Numerator / (float(npts-1))**2.
		  r(ichan,jchan) = Numerator/sqrt((npts*ss(ichan)-s(ichan)*s(ichan))*(npts*ss(jchan)-s(jchan)*s(jchan)))
		 enddo
				
				
		 Write(6,6105) ichan,(cov(ichan,jchan),jchan=1,10)
		 Write(2,6105) ichan,(cov(ichan,jchan),jchan=1,10)
		 Write(6,6105) ichan,(cov(ichan,jchan),jchan=11,20)
		 Write(2,6105) ichan,(cov(ichan,jchan),jchan=11,20)

		enddo

		Do ichan = 1 , 20
		 Write(6,6110) ichan,(r(ichan,jchan),jchan=1,20)
		 Write(2,6110) ichan,(r(ichan,jchan),jchan=1,20)
		EndDo
6110    format(i3,20f7.3)

		stop

310		continue
		write(6,*) 'Read error ',ios
		stop
		end

		Subroutine Correrror(Iseed,Nchan,Gamma,CNoise)

		Integer Iseed	! seed for random number generator, goes both ways
		Integer Nchan	! number of channels
		Real Gamma(Nchan,Nchan) ! the magic gamma multiplier
		Real CNoise(Nchan)		! generated correlated noise

		Real Mean/0.0/, Var/1.0/
		Integer Ichan
		Real Noise(50)  ! make sure this is dimensioned to at least Nchan

		 do ichan = 1 , nchan
		  call normal_random(Iseed,mean,var,noise(Ichan))
		 enddo
		 Call murrv(Nchan,nchan,gamma,nchan,nchan,noise,1,nchan,cnoise)
		
		Return
		End
			SUBROUTINE NORMAL_RANDOM(ISEED,FMEAN,SIGMA,RANDOM)
C
C************************************************************************
C*									*
C*	Module Name:	NORMAL_RANDOM					*
C*									*
C*	Language:	Fortran-77	   Library: [AIMS.UTILITY]	*
C*	Version.Rev:	1.0 26 DEC 85	Programmer: KLEESPIES		*
C*									*
C*	Calling Seq:	CALL NORMAL_RANDOM(ISEED,FMEAN,SIGMA,RANDOM)	*
C*									*
C*	Description:	Computes a random number from a population of	*
C*		random numbers with mean FMEAN and standard deviation	*
C*		SIGMA.  This routine starts with uniformly distributed	*
C*		random numbers on the interval [0,1] and transforms them*
C*		to be randomly distributed with the property 		*
C*		N(FMEAN,SIGMA).  The transformation equation is 	*
C*									*
C*		x = FMEAN + SIGMA*sqrt(-2lny1) * cos(2pi*y2)		*
C*									*
C*		where x is the transformed random variate, and y1 and	*
C*		y2 are two random variates determined by successive 	*
C*		calls to the intrinsic RAN routine.			*
C*									*
C*		A discussion of this transformation may be found in the	*
C*		text "Probability, Statistical Optics, and Data Testing"*
C*		by B.R. Frieden (Springer-Verlag) on page 165.		* 
C*									*
C*	Input Args:	I*4 	ISEED	Seed for random number 		*
C*			R*4	FMEAN	Mean for desired normally 	*
C*					distributed random variate	*
C*			R*4	SIGMA	Standard deviation for desired	*
C*					normally distributed random	*
C*					variate.			*
C*									*
C*	Output Args:	R*4	RANDOM	Random variate from population	*
C*				with mean FMEAN and standard deviation	*
C*				SIGMA.					*
C*									*
C*	Common Blks:	none						*
C*									*
C*	Include:	none						*
C*									*
C*	Externals:	RAN	intrinsic random number generator	*
C*									*
C*	Data Files:	none						*
C*									*
C*	Restrictions:	Note that ISEED is internally reset in this	*
C*			program.  It should be initially set to some	*
C*			large number.  A series of random numbers can	*
C*			be generated without resetting ISEED.		*
C*									*
C*	Error Codes:	none						*
C*									*
C*	Error Messages:	none						*
C*									*
C************************************************************************
C
	INTEGER*4 ISEED
	REAL*4 FMEAN,SIGMA,RANDOM
	REAL*4 Y1,Y2
	REAL*4 TWOPI /6.283185308/
	REAL RAN

c	compute two random variates

  100   CONTINUE

	Y1 = RAN(ISEED)

c	check for zero value (we can't very well go about trying to
c	take the natural log of zero now can we?)

	IF(Y1 .EQ. 0.0) GO TO 100

	Y2 = RAN(ISEED)

	RANDOM = FMEAN + SIGMA*SQRT(-2.0*LOG(Y1))*COS(TWOPI*Y2)

	RETURN
	END