4

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 201x101 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 201x101 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. Can anyone help with this? Thanks in advance!!

Marta Trodahl
  • 41
  • 1
  • 1
    I don't see a programming problem here. Is your question about the algorithm? I would recommend working out the algorithm on paper before programming. In general, it is not a good idea to use lowercase `l` for indexing. – milancurcic Sep 04 '14 at 19:55
  • 2
    For the algorithm scicomp.stackexchange would be better. – Vladimir F Героям слава Sep 04 '14 at 20:15
  • Thanks!! What about the lower case L? Is it because it looks like the number 1, or can it cause an actual problem? I will change it anyways, if it is not a good habit!:) My problem was the algorithm. Stated differently; 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? My code is not summing what I expect, so I am doing something wrong.. – Marta Trodahl Sep 05 '14 at 19:29

1 Answers1

0

First, the number of lattice points that lie within a ring is not a simple function of radius. See the Gauss Circle Problem http://mathworld.wolfram.com/GausssCircleProblem.html

Second, Fortran 90 does have a command similar to the find function in Matlab. It is called where.

Hope that helps.

Norbert S
  • 275
  • 2
  • 14