QMASK Covariance Matrices Fortran Program

        program qmask_demo
        ! Max Tegmark & Angelica de Oliveira-Costa 1999
        ! This code shows how to load a noise covariance matrix file.
        ! Must be compiled with real*8.
        ! For instance, on Solaris, you'd type f77 -r8 qmap_demo.f
        !
        ! The matrices are stored in Solaris fortran binary form.
        ! Since they are symmetric, only the lower triangle is stored
        ! - the subroutines triangular2square and square2triangular
        ! convert between this and a regular square matrix format.
        integer nmax, n, i, j
        parameter(nmax=6500)
        real M(nmax,nmax), sdev1, sdev2, cov, corr
        character*60 fname
        print *,'Number of pixels? (say 2695)'
        read *,n
        if (n.gt.nmax) pause 'DEATH ERROR: nmax too small'
        print *,'Name of covariance matrix file? (say 1Ka12_Sigma.raw)'
        read(*,'(a)'),fname
        call LoadRawVector(n*(n+1)/2,M,fname)
        call triangular2square(nmax,n,M,M)
666     print *,'1st pixel number?'
        read *,i
        print *,'2nd pixel number?'
        read *,j
        if ((i.lt.1).or.(j.lt.1).or.(i.gt.n).or.(j.gt.n)) pause 'Bad pixel number'
        sdev1   = sqrt(M(i,i))
        sdev2   = sqrt(M(j,j))
        cov     = M(i,j)
        corr    = cov/(sdev1*sdev2)
        print *,'RMS noise in 1st pixel.............',sdev1
        print *,'RMS noise in 2nd pixel.............',sdev2
        print *,'Noise covariance ..................',cov
        print *,'Dimensionless noise correlation....',corr
        goto 666
        end
        
        subroutine square2triangular(np,n,A,M)
        ! Converts the n x n symmetric matrix M into my 1-dimensional form.
        ! Because of the FORTRAN number ordering in matrices,
        ! it works even if A and M coincide physically.
        implicit none
        integer np, n, i, j
        real A(np,np), M(n*(n+1)/2)
        do i=1,n
          do j=1,i
            M((i*(i-1))/2 + j) = A(j,i)
          end do
        end do
        return 
        end

        subroutine triangular2square(np,n,M,A)
        ! The inverse of square2triangular.
        ! Unpacks the triangular M into the n x n symmetric matrix A.
        ! Because of the FORTRAN number ordering in matrices,
        ! it works even if A and M coincide physically.
        implicit none
        integer np, n, i, j
        real A(np,np), M(n*(n+1)/2)
        do i=n,1,-1
          do j=i,1,-1
            A(j,i) = M((i*(i-1))/2 + j)
          end do
        end do
        do i=1,n
          do j=1,i-1
            A(i,j) = A(j,i)
          end do
        end do
        return 
        end

        subroutine SaveRawVector(n,a,fname)
        integer n
        real a(n)
        character*60 fname
        open(2,file=fname,form='unformatted')
        print *,'Saving raw vector in ',fname
        write(2) a
        close(2)
        print *,'File size should be ',8*n+8,' bytes.'
        return
        end

        subroutine LoadRawVector(n,a,fname)
        integer n
        real a(n)
        character*60 fname
        open(2,file=fname,form='unformatted',status='old')
        print *,'Loading raw vector from ',fname
        read(2) a
        close(2)
        return
        end

Back To QMASK Product Table