fft2d.F90 3.09 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
MODULE m_fft2d
CONTAINS
  SUBROUTINE fft2d(&
       &                 stars,&
       &                 afft2,bfft2,&
       &                 fg,fgi,fgxy,&
       &                 stride,isn,&
       &                 gfxy )

    !*************************************************************
    !*                                                           *
    !* interface for fg2(star) -- FFT --> gfft (r)     (isn=+1)  *
    !*            or gfft (r)  -- FFT --> fg2(star)    (isn=-1)  *
    !*                                                           *
15
    !* dimension of gfft2 is (3*stars%mx1 x 3*stars%mx2)                     *
16 17 18 19 20 21
    !* afft/bfft contain the real/imaginary part of gfft         *
    !* stars%igfft2(i,1) is the pointer from the G-sphere to stars     *
    !* stars%igfft2(i,2) is the pointer from the G-sphere to fft-grid  *
    !* stars%pgfft2(i)   contains the phases of the G-vectors of sph.  *
    !*                                                           *
    !*************************************************************
22
#include"cpp_double.h"
23 24 25 26 27 28 29
    USE m_cfft
    USE m_types
    IMPLICIT NONE
    TYPE(t_stars),INTENT(IN) :: stars
    INTEGER, INTENT (IN) :: isn,stride
    REAL                 :: fg,fgi

30
    REAL,   INTENT (INOUT):: afft2(0:9*stars%mx1*stars%mx2-1),bfft2(0:9*stars%mx1*stars%mx2-1)
31 32 33 34 35 36 37
    COMPLEX               :: fgxy(stride,stars%ng2-1)
    REAL,OPTIONAL,INTENT(IN) :: gfxy(0:) !factor to calculate the derivates, i.e. g_x

    !... local variables

    INTEGER i,ifftd2
    REAL  scale
38
    COMPLEX fg2(stars%ng2)
39

40
    ifftd2=9*stars%mx1*stars%mx2
41 42 43 44 45 46
    !
    IF (isn.GT.0) THEN
       !
       !  ---> put stars onto the fft-grid 
       !
       fg2(1) = CMPLX(fg,fgi)
47 48
       CALL CPP_BLAS_ccopy(stars%ng2-1,fgxy,stride,fg2(2),1)
       !fg2(2:)=fgxy(1,:)
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
       afft2=0.0
       bfft2=0.0
       IF (PRESENT(gfxy)) THEN
          DO i=0,stars%kimax2
             afft2(stars%igfft2(i,2))=REAL(fg2(stars%igfft2(i,1))*stars%pgfft2(i))*gfxy(i)
             bfft2(stars%igfft2(i,2))=AIMAG(fg2(stars%igfft2(i,1))*stars%pgfft2(i))*gfxy(i)
          ENDDO
       ELSE 
          DO i=0,stars%kimax2
             afft2(stars%igfft2(i,2))=REAL(fg2(stars%igfft2(i,1))*stars%pgfft2(i))
             bfft2(stars%igfft2(i,2))=AIMAG(fg2(stars%igfft2(i,1))*stars%pgfft2(i))
          ENDDO
       ENDIF
    ENDIF

    !---> now do the fft (isn=+1 : G -> r ; isn=-1 : r -> G)

66 67
    CALL cfft(afft2,bfft2,ifftd2,3*stars%mx1,3*stars%mx1,isn)
    CALL cfft(afft2,bfft2,ifftd2,3*stars%mx2,ifftd2,isn)
68 69 70 71 72 73 74 75 76 77

    IF (isn.LT.0) THEN
       !
       !  ---> collect stars from the fft-grid
       !
       DO i=1,stars%ng2
          fg2(i) = CMPLX(0.0,0.0)
       ENDDO
       scale=1.0/ifftd2
       DO i=0,stars%kimax2
78
          fg2(stars%igfft2(i,1))=fg2(stars%igfft2(i,1))+ CONJG( stars%pgfft2(i) ) * &
79 80 81 82 83 84 85 86 87 88 89
               &                 CMPLX(afft2(stars%igfft2(i,2)),bfft2(stars%igfft2(i,2)))
       ENDDO
       fg=scale*REAL(fg2(1))/stars%nstr2(1)
       fgi=scale*AIMAG(fg2(1))/stars%nstr2(1)
       DO i=2,stars%ng2
          fgxy(1,i-1)=scale*fg2(i)/stars%nstr2(i)
       ENDDO
    ENDIF

  END SUBROUTINE fft2d
END MODULE m_fft2d