Azimuthal average in Fortran? Find indexes in Fortran?

Question

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

!Radial distance

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? )

Answer 1

I haven't spent much time looking at your code, but I would approach this in a slightly different way. Rather than asking which wave numbers correspond to a given radius range, I would divide the space into the (radial) bins you want and then loop over the wave vectors asking which bin contains each wave vector (e.g. using integer division). By summing the values in each bin and counting how many are in each bin you can easily compute the averages. This should be more efficient too.

Answer 2

To find indices in Fortran arrays: create your sellecting mask, and use MAXLOC or MINLOC