I am working on an eigenvalue problem in fortran. I have used Lapack to solve the problem and get the eigenvalues and eigenvectors. This is done for $201\times101$ wavenumbers, only half the wavespace due to symmetry, and for each grid-point in a large domain(in the ocean). I am searching for the maximum eigenvalue for each gridpoint, and I would like to not just pick the absolute maximum in the $201\times101$ matrix of eigenvalues, but perform an azimuthal average in wavespace - and then pick the maximum average. I am struggeling with seeing how to do this.
At first I coded it like this:
! Wavenumber domain
dx=4000. pi = 4.*atan(1.) DO m=1,ktot kx(m) = -(2.*pi)/(dx) + ((m-1)*2.*pi)/(100.*dx) END DO DO l=1,ltot ly(l)= ((l-1)*2.*pi)/(100.*dx) END DO
DO m=1,ktot DO l=1,ltot raddist(m,l)=sqrt(kx(m)**2+ly(l)**2) END DO END DO
! Azimuthal average (omegai are the eigenvalues I have found, a
ktot*ltot large matrix)
DO i=1,ltot-1 ind=(raddist(:,i).GE.ly(i).AND.raddist(:,i).LT.ly(i+1)) length=count(ind) WHERE (ind) average_omegai = sum(omegai)/length END DO
But it seems I am summing a horizontal chunk for all wavenumbers in k-direction, between
ly(i) and ly(i+1). I rather need to make a semicircle in wavespace to sum all omegai-values in between it. Can anyone help with this?
Thanks so much in advance!!
(BTW: I also tried
DO m=1,ktot-1 DO l=1,ltot-1 ind=(raddist(m,l).GE.ly(l).AND.raddist(m,l).LT.ly(l+1)) length=count(ind) WHERE (ind) average_omegai = sum(omegai)/length END DO END DO
but I right now I am blind to what I am doing wrong! My problem stated differently would be: is there another way in Fortran to find indexes in a matrix which fulfill some criteria? The find-function in Matlab would do the work, but I have not found an equivalent intrinsic function in Fortran? )