  	Program AMSUA_Corr_Noise

!	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 /15/
		
		Integer Iseed /72338/ ! just a number I made up.
		Real cnoise(15,1000000)
		real ggamma(15,15)

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

		Real Cov(15,15)
		Real r(15,15)

		Integer Npts
		Real*8 elapsed_time

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

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

	    
!		write(6,*) '11 or 12'
!		read(5,*) instr
!		if(instr .eq. 11) then
!		 open(3,file='..\h11ggamma.dat',status='old',form='formatted')    
!		else
!		 open(3,file='..\h12ggamma.dat',status='old',form='formatted')
!		endif
        write(6,*) '15 or 16'
		read(5,*) instr
		if(instr .eq. 15) then
		 open(3,file='..\AMSUA_SpaceStats\N15_AMSUAgammanew.dat',status='old',form='formatted')    
		else
		 open(3,file='..\AMSUA_SpaceStats\N16_AMSUAgammanew.dat',status='old',form='formatted')    
		endif


		do k = 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,15
		   s(ichan) = s(ichan) + cnoise(ichan,k)
          ss(ichan) =ss(ichan) + cnoise(ichan,k)*cnoise(ichan,k)
		  do jchan = 1 , 15
		   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 , 15
		 do jchan = 1 , 15
		  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,15)
		 Write(2,6105) ichan,(cov(ichan,jchan),jchan=11,15)

		enddo

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

		stop

310		continue
		write(6,*) 'Read error ',ios
		stop
End Program AMSUA_Corr_Noise

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

	INTEGER*4 ISEED
	REAL*4 FMEAN,SIGMA,RANDOM
	REAL*4 Y1,Y2
	REAL*4 TWOPI /6.283185308/
	REAL RAN

!	compute two random variates

  100   CONTINUE

	Y1 = RAN(ISEED)

!	check for zero value (we can't very well go about trying to
!	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

