text:
! Fast Fourier/Cosine/Sine Transform
!     dimension   :one
!     data length :power of 2
!     decimation  :frequency
!     radix       :8, 4, 2
!     data        :inplace
!     table       :use
! subroutines
!     cdft: Complex Discrete Fourier Transform
!     rdft: Real Discrete Fourier Transform
!     ddct: Discrete Cosine Transform
!     ddst: Discrete Sine Transform
!     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
!     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
!
!
! -------- Complex DFT (Discrete Fourier Transform) --------
!     [definition]
!         <case1>
!             X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n
!         <case2>
!             X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n
!         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call cdft(2*n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call cdft(2*n, -1, a, ip, w)
!     [parameters]
!         2*n          :data length (integer)
!                       n >= 1, n = power of 2
!         a(0:2*n-1)   :input/output data (real*8)
!                       input data
!                           a(2*j) = Re(x(j)),
!                           a(2*j+1) = Im(x(j)), 0<=j<n
!                       output data
!                           a(2*k) = Re(X(k)),
!                           a(2*k+1) = Im(X(k)), 0<=k<n
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n/2-1)   :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call cdft(2*n, -1, a, ip, w)
!         is
!             call cdft(2*n, 1, a, ip, w)
!             do j = 0, 2 * n - 1
!                 a(j) = a(j) / n
!             end do
!         .
!
!
! -------- Real DFT / Inverse of Real DFT --------
!     [definition]
!         <case1> RDFT
!             R(k) = sum_j=0^n-1 a(j)*cos(2*pi*j*k/n), 0<=k<=n/2
!             I(k) = sum_j=0^n-1 a(j)*sin(2*pi*j*k/n), 0<k<n/2
!         <case2> IRDFT (excluding scale)
!             a(k) = (R(0) + R(n/2)*cos(pi*k))/2 +
!                    sum_j=1^n/2-1 R(j)*cos(2*pi*j*k/n) +
!                    sum_j=1^n/2-1 I(j)*sin(2*pi*j*k/n), 0<=k<n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call rdft(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call rdft(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       <case1>
!                           output data
!                               a(2*k) = R(k), 0<=k<n/2
!                               a(2*k+1) = I(k), 0<k<n/2
!                               a(1) = R(n/2)
!                       <case2>
!                           input data
!                               a(2*j) = R(j), 0<=j<n/2
!                               a(2*j+1) = I(j), 0<j<n/2
!                               a(1) = R(n/2)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n/2-1)   :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call rdft(n, 1, a, ip, w)
!         is
!             call rdft(n, -1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
!     [definition]
!         <case1> IDCT (excluding scale)
!             C(k) = sum_j=0^n-1 a(j)*cos(pi*j*(k+1/2)/n), 0<=k<n
!         <case2> DCT
!             C(k) = sum_j=0^n-1 a(j)*cos(pi*(j+1/2)*k/n), 0<=k<n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddct(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddct(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       output data
!                           a(k) = C(k), 0<=k<n
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/4-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call ddct(n, -1, a, ip, w)
!         is
!             a(0) = a(0) / 2
!             call ddct(n, 1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- DST (Discrete Sine Transform) / Inverse of DST --------
!     [definition]
!         <case1> IDST (excluding scale)
!             S(k) = sum_j=1^n A(j)*sin(pi*j*(k+1/2)/n), 0<=k<n
!         <case2> DST
!             S(k) = sum_j=0^n-1 a(j)*sin(pi*(j+1/2)*k/n), 0<k<=n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddst(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddst(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       <case1>
!                           input data
!                               a(j) = A(j), 0<j<n
!                               a(0) = A(n)
!                           output data
!                               a(k) = S(k), 0<=k<n
!                       <case2>
!                           output data
!                               a(k) = S(k), 0<k<n
!                               a(0) = S(n)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/4-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call ddst(n, -1, a, ip, w)
!         is
!             a(0) = a(0) / 2
!             call ddst(n, 1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
!     [definition]
!         C(k) = sum_j=0^n a(j)*cos(pi*j*k/n), 0<=k<=n
!     [usage]
!         ip(0) = 0  ! first time only
!         call dfct(n, a, t, ip, w)
!     [parameters]
!         n            :data length - 1 (integer)
!                       n >= 2, n = power of 2
!         a(0:n)       :input/output data (real*8)
!                       output data
!                           a(k) = C(k), 0<=k<=n
!         t(0:n/2)     :work area (real*8)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/4)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/4+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/8-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             a(0) = a(0) / 2
!             a(n) = a(n) / 2
!             call dfct(n, a, t, ip, w)
!         is
!             a(0) = a(0) / 2
!             a(n) = a(n) / 2
!             call dfct(n, a, t, ip, w)
!             do j = 0, n
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
!     [definition]
!         S(k) = sum_j=1^n-1 a(j)*sin(pi*j*k/n), 0<k<n
!     [usage]
!         ip(0) = 0  ! first time only
!         call dfst(n, a, t, ip, w)
!     [parameters]
!         n            :data length + 1 (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       output data
!                           a(k) = S(k), 0<k<n
!                       (a(0) is used for work area)
!         t(0:n/2-1)   :work area (real*8)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/4)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/4+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/8-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call dfst(n, a, t, ip, w)
!         is
!             call dfst(n, a, t, ip, w)
!             do j = 1, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! Appendix :
!     The cos/sin table is recalculated when the larger table required.
!     w() and ip() are compatible with all routines.
!
!
      subroutine cdft(n, isgn, a, ip, w)
      integer n, isgn, ip(0 : *)
      real*8 a(0 : n - 1), w(0 : *)
      if (n .gt. 4 * ip(0)) then
          call makewt(n / 4, ip, w)
      end if
      if (n .gt. 4) then
          if (isgn .ge. 0) then
              call bitrv2(n, ip(2), a)
              call cftfsub(n, a, w)
          else
              call bitrv2conj(n, ip(2), a)
              call cftbsub(n, a, w)
          end if
      else if (n .eq. 4) then
          call cftfsub(n, a, w)
      end if
      end
!
      subroutine rdft(n, isgn, a, ip, w)
      integer n, isgn, ip(0 : *), nw, nc
      real*8 a(0 : n - 1), w(0 : *), xi
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 4 * nc) then
          nc = n / 4
          call makect(nc, ip, w(nw))
      end if
      if (isgn .ge. 0) then
          if (n .gt. 4) then
              call bitrv2(n, ip(2), a)
              call cftfsub(n, a, w)
              call rftfsub(n, a, nc, w(nw))
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
          xi = a(0) - a(1)
          a(0) = a(0) + a(1)
          a(1) = xi
      else
          a(1) = 0.5d0 * (a(0) - a(1))
          a(0) = a(0) - a(1)
          if (n .gt. 4) then
              call rftbsub(n, a, nc, w(nw))
              call bitrv2(n, ip(2), a)
              call cftbsub(n, a, w)
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
      end if
      end
!
      subroutine ddct(n, isgn, a, ip, w)
      integer n, isgn, ip(0 : *), j, nw, nc
      real*8 a(0 : n - 1), w(0 : *), xr
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
          xr = a(n - 1)
          do j = n - 2, 2, -2
              a(j + 1) = a(j) - a(j - 1)
              a(j) = a(j) + a(j - 1)
          end do
          a(1) = a(0) - xr
          a(0) = a(0) + xr
          if (n .gt. 4) then
              call rftbsub(n, a, nc, w(nw))
              call bitrv2(n, ip(2), a)
              call cftbsub(n, a, w)
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
      end if
      call dctsub(n, a, nc, w(nw))
      if (isgn .ge. 0) then
          if (n .gt. 4) then
              call bitrv2(n, ip(2), a)
              call cftfsub(n, a, w)
              call rftfsub(n, a, nc, w(nw))
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
          xr = a(0) - a(1)
          a(0) = a(0) + a(1)
          do j = 2, n - 2, 2
              a(j - 1) = a(j) - a(j + 1)
              a(j) = a(j) + a(j + 1)
          end do
          a(n - 1) = xr
      end if
      end
!
      subroutine ddst(n, isgn, a, ip, w)
      integer n, isgn, ip(0 : *), j, nw, nc
      real*8 a(0 : n - 1), w(0 : *), xr
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
          xr = a(n - 1)
          do j = n - 2, 2, -2
              a(j + 1) = -a(j) - a(j - 1)
              a(j) = a(j) - a(j - 1)
          end do
          a(1) = a(0) + xr
          a(0) = a(0) - xr
          if (n .gt. 4) then
              call rftbsub(n, a, nc, w(nw))
              call bitrv2(n, ip(2), a)
              call cftbsub(n, a, w)
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
      end if
      call dstsub(n, a, nc, w(nw))
      if (isgn .ge. 0) then
          if (n .gt. 4) then
              call bitrv2(n, ip(2), a)
              call cftfsub(n, a, w)
              call rftfsub(n, a, nc, w(nw))
          else if (n .eq. 4) then
              call cftfsub(n, a, w)
          end if
          xr = a(0) - a(1)
          a(0) = a(0) + a(1)
          do j = 2, n - 2, 2
              a(j - 1) = -a(j) - a(j + 1)
              a(j) = a(j) - a(j + 1)
          end do
          a(n - 1) = -xr
      end if
      end
!
      subroutine dfct(n, a, t, ip, w)
      integer n, ip(0 : *), j, k, l, m, mh, nw, nc
      real*8 a(0 : n), t(0 : n / 2), w(0 : *), xr, xi, yr, yi
      nw = ip(0)
      if (n .gt. 8 * nw) then
          nw = n / 8
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 2 * nc) then
          nc = n / 2
          call makect(nc, ip, w(nw))
      end if
      m = n / 2
      yi = a(m)
      xi = a(0) + a(n)
      a(0) = a(0) - a(n)
      t(0) = xi - yi
      t(m) = xi + yi
      if (n .gt. 2) then
          mh = m / 2
          do j = 1, mh - 1
              k = m - j
              xr = a(j) - a(n - j)
              xi = a(j) + a(n - j)
              yr = a(k) - a(n - k)
              yi = a(k) + a(n - k)
              a(j) = xr
              a(k) = yr
              t(j) = xi - yi
              t(k) = xi + yi
          end do
          t(mh) = a(mh) + a(n - mh)
          a(mh) = a(mh) - a(n - mh)
          call dctsub(m, a, nc, w(nw))
          if (m .gt. 4) then
              call bitrv2(m, ip(2), a)
              call cftfsub(m, a, w)
              call rftfsub(m, a, nc, w(nw))
          else if (m .eq. 4) then
              call cftfsub(m, a, w)
          end if
          a(n - 1) = a(0) - a(1)
          a(1) = a(0) + a(1)
          do j = m - 2, 2, -2
              a(2 * j + 1) = a(j) + a(j + 1)
              a(2 * j - 1) = a(j) - a(j + 1)
          end do
          l = 2
          m = mh
          do while (m .ge. 2)
              call dctsub(m, t, nc, w(nw))
              if (m .gt. 4) then
                  call bitrv2(m, ip(2), t)
                  call cftfsub(m, t, w)
                  call rftfsub(m, t, nc, w(nw))
              else if (m .eq. 4) then
                  call cftfsub(m, t, w)
              end if
              a(n - l) = t(0) - t(1)
              a(l) = t(0) + t(1)
              k = 0
              do j = 2, m - 2, 2
                  k = k + 4 * l
                  a(k - l) = t(j) - t(j + 1)
                  a(k + l) = t(j) + t(j + 1)
              end do
              l = 2 * l
              mh = m / 2
              do j = 0, mh - 1
                  k = m - j
                  t(j) = t(m + k) - t(m + j)
                  t(k) = t(m + k) + t(m + j)
              end do
              t(mh) = t(m + mh)
              m = mh
          end do
          a(l) = t(0)
          a(n) = t(2) - t(1)
          a(0) = t(2) + t(1)
      else
          a(1) = a(0)
          a(2) = t(0)
          a(0) = t(1)
      end if
      end
!
      subroutine dfst(n, a, t, ip, w)
      integer n, ip(0 : *), j, k, l, m, mh, nw, nc
      real*8 a(0 : n - 1), t(0 : n / 2 - 1), w(0 : *), xr, xi, yr, yi
      nw = ip(0)
      if (n .gt. 8 * nw) then
          nw = n / 8
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 2 * nc) then
          nc = n / 2
          call makect(nc, ip, w(nw))
      end if
      if (n .gt. 2) then
          m = n / 2
          mh = m / 2
          do j = 1, mh - 1
              k = m - j
              xr = a(j) + a(n - j)
              xi = a(j) - a(n - j)
              yr = a(k) + a(n - k)
              yi = a(k) - a(n - k)
              a(j) = xr
              a(k) = yr
              t(j) = xi + yi
              t(k) = xi - yi
          end do
          t(0) = a(mh) - a(n - mh)
          a(mh) = a(mh) + a(n - mh)
          a(0) = a(m)
          call dstsub(m, a, nc, w(nw))
          if (m .gt. 4) then
              call bitrv2(m, ip(2), a)
              call cftfsub(m, a, w)
              call rftfsub(m, a, nc, w(nw))
          else if (m .eq. 4) then
              call cftfsub(m, a, w)
          end if
          a(n - 1) = a(1) - a(0)
          a(1) = a(0) + a(1)
          do j = m - 2, 2, -2
              a(2 * j + 1) = a(j) - a(j + 1)
              a(2 * j - 1) = -a(j) - a(j + 1)
          end do
          l = 2
          m = mh
          do while (m .ge. 2)
              call dstsub(m, t, nc, w(nw))
              if (m .gt. 4) then
                  call bitrv2(m, ip(2), t)
                  call cftfsub(m, t, w)
                  call rftfsub(m, t, nc, w(nw))
              else if (m .eq. 4) then
                  call cftfsub(m, t, w)
              end if
              a(n - l) = t(1) - t(0)
              a(l) = t(0) + t(1)
              k = 0
              do j = 2, m - 2, 2
                  k = k + 4 * l
                  a(k - l) = -t(j) - t(j + 1)
                  a(k + l) = t(j) - t(j + 1)
              end do
              l = 2 * l
              mh = m / 2
              do j = 1, mh - 1
                  k = m - j
                  t(j) = t(m + k) + t(m + j)
                  t(k) = t(m + k) - t(m + j)
              end do
              t(0) = t(m + mh)
              m = mh
          end do
          a(l) = t(0)
      end if
      a(0) = 0
      end
!
! -------- initializing routines --------
!
      subroutine makewt(nw, ip, w)
      integer nw, ip(0 : *), j, nwh
      real*8 w(0 : nw - 1), delta, x, y
      ip(0) = nw
      ip(1) = 1
      if (nw .gt. 2) then
          nwh = nw / 2
          delta = atan(1.0d0) / nwh
          w(0) = 1
          w(1) = 0
          w(nwh) = cos(delta * nwh)
          w(nwh + 1) = w(nwh)
          if (nwh .gt. 2) then
              do j = 2, nwh - 2, 2
                  x = cos(delta * j)
                  y = sin(delta * j)
                  w(j) = x
                  w(j + 1) = y
                  w(nw - j) = y
                  w(nw - j + 1) = x
              end do
              do j = nwh - 2, 2, -2
                  x = w(2 * j)
                  y = w(2 * j + 1)
                  w(nwh + j) = x
                  w(nwh + j + 1) = y
              end do
              call bitrv2(nw, ip(2), w)
          end if
      end if
      end
!
      subroutine makect(nc, ip, c)
      integer nc, ip(0 : *), j, nch
      real*8 c(0 : nc - 1), delta
      ip(1) = nc
      if (nc .gt. 1) then
          nch = nc / 2
          delta = atan(1.0d0) / nch
          c(0) = cos(delta * nch)
          c(nch) = 0.5d0 * c(0)
          do j = 1, nch - 1
              c(j) = 0.5d0 * cos(delta * j)
              c(nc - j) = 0.5d0 * sin(delta * j)
          end do
      end if
      end
!
! -------- child routines --------
!
      subroutine bitrv2(n, ip, a)
      integer n, ip(0 : *), j, j1, k, k1, l, m, m2
      real*8 a(0 : n - 1), xr, xi, yr, yi
      ip(0) = 0
      l = n
      m = 1
      do while (8 * m .lt. l)
          l = l / 2
          do j = 0, m - 1
              ip(m + j) = ip(j) + l
          end do
          m = m * 2
      end do
      m2 = 2 * m
      if (8 * m .eq. l) then
          do k = 0, m - 1
              do j = 0, k - 1
                  j1 = 2 * j + ip(k)
                  k1 = 2 * k + ip(j)
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + 2 * m2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 - m2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + 2 * m2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
              end do
              j1 = 2 * k + m2 + ip(k)
              k1 = j1 + m2
              xr = a(j1)
              xi = a(j1 + 1)
              yr = a(k1)
              yi = a(k1 + 1)
              a(j1) = yr
              a(j1 + 1) = yi
              a(k1) = xr
              a(k1 + 1) = xi
          end do
      else
          do k = 1, m - 1
              do j = 0, k - 1
                  j1 = 2 * j + ip(k)
                  k1 = 2 * k + ip(j)
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + m2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
              end do
          end do
      end if
      end
!
      subroutine bitrv2conj(n, ip, a)
      integer n, ip(0 : *), j, j1, k, k1, l, m, m2
      real*8 a(0 : n - 1), xr, xi, yr, yi
      ip(0) = 0
      l = n
      m = 1
      do while (8 * m .lt. l)
          l = l / 2
          do j = 0, m - 1
              ip(m + j) = ip(j) + l
          end do
          m = m * 2
      end do
      m2 = 2 * m
      if (8 * m .eq. l) then
          do k = 0, m - 1
              do j = 0, k - 1
                  j1 = 2 * j + ip(k)
                  k1 = 2 * k + ip(j)
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + 2 * m2
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 - m2
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + 2 * m2
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
              end do
              k1 = 2 * k + ip(k)
              a(k1 + 1) = -a(k1 + 1)
              j1 = k1 + m2
              k1 = j1 + m2
              xr = a(j1)
              xi = -a(j1 + 1)
              yr = a(k1)
              yi = -a(k1 + 1)
              a(j1) = yr
              a(j1 + 1) = yi
              a(k1) = xr
              a(k1 + 1) = xi
              k1 = k1 + m2
              a(k1 + 1) = -a(k1 + 1)
          end do
      else
          a(1) = -a(1)
          a(m2 + 1) = -a(m2 + 1)
          do k = 1, m - 1
              do j = 0, k - 1
                  j1 = 2 * j + ip(k)
                  k1 = 2 * k + ip(j)
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + m2
                  k1 = k1 + m2
                  xr = a(j1)
                  xi = -a(j1 + 1)
                  yr = a(k1)
                  yi = -a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
              end do
              k1 = 2 * k + ip(k)
              a(k1 + 1) = -a(k1 + 1)
              a(k1 + m2 + 1) = -a(k1 + m2 + 1)
          end do
      end if
      end
!
      subroutine cftfsub(n, a, w)
      integer n, j, j1, j2, j3, l
      real*8 a(0 : n - 1), w(0 : *)
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      l = 2
      if (n .ge. 16) then
          call cft1st(n, a, w)
          l = 16
          do while (8 * l .le. n)
              call cftmdl(n, l, a, w)
              l = 8 * l
          end do
      end if
      if (2 * l .lt. n) then
          do j = 0, l - 2, 2
              j1 = j + l
              j2 = j1 + l
              j3 = j2 + l
              x0r = a(j) + a(j1)
              x0i = a(j + 1) + a(j1 + 1)
              x1r = a(j) - a(j1)
              x1i = a(j + 1) - a(j1 + 1)
              x2r = a(j2) + a(j3)
              x2i = a(j2 + 1) + a(j3 + 1)
              x3r = a(j2) - a(j3)
              x3i = a(j2 + 1) - a(j3 + 1)
              a(j) = x0r + x2r
              a(j + 1) = x0i + x2i
              a(j2) = x0r - x2r
              a(j2 + 1) = x0i - x2i
              a(j1) = x1r - x3i
              a(j1 + 1) = x1i + x3r
              a(j3) = x1r + x3i
              a(j3 + 1) = x1i - x3r
          end do
      else if (2 * l .eq. n) then
          do j = 0, l - 2, 2
              j1 = j + l
              x0r = a(j) - a(j1)
              x0i = a(j + 1) - a(j1 + 1)
              a(j) = a(j) + a(j1)
              a(j + 1) = a(j + 1) + a(j1 + 1)
              a(j1) = x0r
              a(j1 + 1) = x0i
          end do
      end if
      end
!
      subroutine cftbsub(n, a, w)
      integer n, j, j1, j2, j3, j4, j5, j6, j7, l
      real*8 a(0 : n - 1), w(0 : *)
      real*8 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      l = 2
      if (n .gt. 16) then
          call cft1st(n, a, w)
          l = 16
          do while (8 * l .lt. n)
              call cftmdl(n, l, a, w)
              l = 8 * l
          end do
      end if
      if (4 * l .lt. n) then
          wn4r = w(2)
          do j = 0, l - 2, 2
              j1 = j + l
              j2 = j1 + l
              j3 = j2 + l
              j4 = j3 + l
              j5 = j4 + l
              j6 = j5 + l
              j7 = j6 + l
              x0r = a(j) + a(j1)
              x0i = -a(j + 1) - a(j1 + 1)
              x1r = a(j) - a(j1)
              x1i = -a(j + 1) + a(j1 + 1)
              x2r = a(j2) + a(j3)
              x2i = a(j2 + 1) + a(j3 + 1)
              x3r = a(j2) - a(j3)
              x3i = a(j2 + 1) - a(j3 + 1)
              y0r = x0r + x2r
              y0i = x0i - x2i
              y2r = x0r - x2r
              y2i = x0i + x2i
              y1r = x1r - x3i
              y1i = x1i - x3r
              y3r = x1r + x3i
              y3i = x1i + x3r
              x0r = a(j4) + a(j5)
              x0i = a(j4 + 1) + a(j5 + 1)
              x1r = a(j4) - a(j5)
              x1i = a(j4 + 1) - a(j5 + 1)
              x2r = a(j6) + a(j7)
              x2i = a(j6 + 1) + a(j7 + 1)
              x3r = a(j6) - a(j7)
              x3i = a(j6 + 1) - a(j7 + 1)
              y4r = x0r + x2r
              y4i = x0i + x2i
              y6r = x0r - x2r
              y6i = x0i - x2i
              x0r = x1r - x3i
              x0i = x1i + x3r
              x2r = x1r + x3i
              x2i = x1i - x3r
              y5r = wn4r * (x0r - x0i)
              y5i = wn4r * (x0r + x0i)
              y7r = wn4r * (x2r - x2i)
              y7i = wn4r * (x2r + x2i)
              a(j1) = y1r + y5r
              a(j1 + 1) = y1i - y5i
              a(j5) = y1r - y5r
              a(j5 + 1) = y1i + y5i
              a(j3) = y3r - y7i
              a(j3 + 1) = y3i - y7r
              a(j7) = y3r + y7i
              a(j7 + 1) = y3i + y7r
              a(j) = y0r + y4r
              a(j + 1) = y0i - y4i
              a(j4) = y0r - y4r
              a(j4 + 1) = y0i + y4i
              a(j2) = y2r - y6i
              a(j2 + 1) = y2i - y6r
              a(j6) = y2r + y6i
              a(j6 + 1) = y2i + y6r
          end do
      else if (4 * l .eq. n) then
          do j = 0, l - 2, 2
              j1 = j + l
              j2 = j1 + l
              j3 = j2 + l
              x0r = a(j) + a(j1)
              x0i = -a(j + 1) - a(j1 + 1)
              x1r = a(j) - a(j1)
              x1i = -a(j + 1) + a(j1 + 1)
              x2r = a(j2) + a(j3)
              x2i = a(j2 + 1) + a(j3 + 1)
              x3r = a(j2) - a(j3)
              x3i = a(j2 + 1) - a(j3 + 1)
              a(j) = x0r + x2r
              a(j + 1) = x0i - x2i
              a(j2) = x0r - x2r
              a(j2 + 1) = x0i + x2i
              a(j1) = x1r - x3i
              a(j1 + 1) = x1i - x3r
              a(j3) = x1r + x3i
              a(j3 + 1) = x1i + x3r
          end do
      else
          do j = 0, l - 2, 2
              j1 = j + l
              x0r = a(j) - a(j1)
              x0i = -a(j + 1) + a(j1 + 1)
              a(j) = a(j) + a(j1)
              a(j + 1) = -a(j + 1) - a(j1 + 1)
              a(j1) = x0r
              a(j1 + 1) = x0i
          end do
      end if
      end
!
      subroutine cft1st(n, a, w)
      integer n, j, k1
      real*8 a(0 : n - 1), w(0 : *)
      real*8 wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i
      real*8 wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      wn4r = w(2)
      x0r = a(0) + a(2)
      x0i = a(1) + a(3)
      x1r = a(0) - a(2)
      x1i = a(1) - a(3)
      x2r = a(4) + a(6)
      x2i = a(5) + a(7)
      x3r = a(4) - a(6)
      x3i = a(5) - a(7)
      y0r = x0r + x2r
      y0i = x0i + x2i
      y2r = x0r - x2r
      y2i = x0i - x2i
      y1r = x1r - x3i
      y1i = x1i + x3r
      y3r = x1r + x3i
      y3i = x1i - x3r
      x0r = a(8) + a(10)
      x0i = a(9) + a(11)
      x1r = a(8) - a(10)
      x1i = a(9) - a(11)
      x2r = a(12) + a(14)
      x2i = a(13) + a(15)
      x3r = a(12) - a(14)
      x3i = a(13) - a(15)
      y4r = x0r + x2r
      y4i = x0i + x2i
      y6r = x0r - x2r
      y6i = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      x2r = x1r + x3i
      x2i = x1i - x3r
      y5r = wn4r * (x0r - x0i)
      y5i = wn4r * (x0r + x0i)
      y7r = wn4r * (x2r - x2i)
      y7i = wn4r * (x2r + x2i)
      a(2) = y1r + y5r
      a(3) = y1i + y5i
      a(10) = y1r - y5r
      a(11) = y1i - y5i
      a(6) = y3r - y7i
      a(7) = y3i + y7r
      a(14) = y3r + y7i
      a(15) = y3i - y7r
      a(0) = y0r + y4r
      a(1) = y0i + y4i
      a(8) = y0r - y4r
      a(9) = y0i - y4i
      a(4) = y2r - y6i
      a(5) = y2i + y6r
      a(12) = y2r + y6i
      a(13) = y2i - y6r
      if (n .gt. 16) then
          wk1r = w(4)
          wk1i = w(5)
          x0r = a(16) + a(18)
          x0i = a(17) + a(19)
          x1r = a(16) - a(18)
          x1i = a(17) - a(19)
          x2r = a(20) + a(22)
          x2i = a(21) + a(23)
          x3r = a(20) - a(22)
          x3i = a(21) - a(23)
          y0r = x0r + x2r
          y0i = x0i + x2i
          y2r = x0r - x2r
          y2i = x0i - x2i
          y1r = x1r - x3i
          y1i = x1i + x3r
          y3r = x1r + x3i
          y3i = x1i - x3r
          x0r = a(24) + a(26)
          x0i = a(25) + a(27)
          x1r = a(24) - a(26)
          x1i = a(25) - a(27)
          x2r = a(28) + a(30)
          x2i = a(29) + a(31)
          x3r = a(28) - a(30)
          x3i = a(29) - a(31)
          y4r = x0r + x2r
          y4i = x0i + x2i
          y6r = x0r - x2r
          y6i = x0i - x2i
          x0r = x1r - x3i
          x0i = x1i + x3r
          x2r = x1r + x3i
          x2i = x3r - x1i
          y5r = wk1i * x0r - wk1r * x0i
          y5i = wk1i * x0i + wk1r * x0r
          y7r = wk1r * x2r + wk1i * x2i
          y7i = wk1r * x2i - wk1i * x2r
          x0r = wk1r * y1r - wk1i * y1i
          x0i = wk1r * y1i + wk1i * y1r
          a(18) = x0r + y5r
          a(19) = x0i + y5i
          a(26) = y5i - x0i
          a(27) = x0r - y5r
          x0r = wk1i * y3r - wk1r * y3i
          x0i = wk1i * y3i + wk1r * y3r
          a(22) = x0r - y7r
          a(23) = x0i + y7i
          a(30) = y7i - x0i
          a(31) = x0r + y7r
          a(16) = y0r + y4r
          a(17) = y0i + y4i
          a(24) = y4i - y0i
          a(25) = y0r - y4r
          x0r = y2r - y6i
          x0i = y2i + y6r
          a(20) = wn4r * (x0r - x0i)
          a(21) = wn4r * (x0i + x0r)
          x0r = y6r - y2i
          x0i = y2r + y6i
          a(28) = wn4r * (x0r - x0i)
          a(29) = wn4r * (x0i + x0r)
          k1 = 4
          do j = 32, n - 16, 16
              k1 = k1 + 4
              wk1r = w(k1)
              wk1i = w(k1 + 1)
              wk2r = w(k1 + 2)
              wk2i = w(k1 + 3)
              wtmp = 2 * wk2i
              wk3r = wk1r - wtmp * wk1i
              wk3i = wtmp * wk1r - wk1i
              wk4r = 1 - wtmp * wk2i
              wk4i = wtmp * wk2r
              wtmp = 2 * wk4i
              wk5r = wk3r - wtmp * wk1i
              wk5i = wtmp * wk1r - wk3i
              wk6r = wk2r - wtmp * wk2i
              wk6i = wtmp * wk2r - wk2i
              wk7r = wk1r - wtmp * wk3i
              wk7i = wtmp * wk3r - wk1i
              x0r = a(j) + a(j + 2)
              x0i = a(j + 1) + a(j + 3)
              x1r = a(j) - a(j + 2)
              x1i = a(j + 1) - a(j + 3)
              x2r = a(j + 4) + a(j + 6)
              x2i = a(j + 5) + a(j + 7)
              x3r = a(j + 4) - a(j + 6)
              x3i = a(j + 5) - a(j + 7)
              y0r = x0r + x2r
              y0i = x0i + x2i
              y2r = x0r - x2r
              y2i = x0i - x2i
              y1r = x1r - x3i
              y1i = x1i + x3r
              y3r = x1r + x3i
              y3i = x1i - x3r
              x0r = a(j + 8) + a(j + 10)
              x0i = a(j + 9) + a(j + 11)
              x1r = a(j + 8) - a(j + 10)
              x1i = a(j + 9) - a(j + 11)
              x2r = a(j + 12) + a(j + 14)
              x2i = a(j + 13) + a(j + 15)
              x3r = a(j + 12) - a(j + 14)
              x3i = a(j + 13) - a(j + 15)
              y4r = x0r + x2r
              y4i = x0i + x2i
              y6r = x0r - x2r
              y6i = x0i - x2i
              x0r = x1r - x3i
              x0i = x1i + x3r
              x2r = x1r + x3i
              x2i = x1i - x3r
              y5r = wn4r * (x0r - x0i)
              y5i = wn4r * (x0r + x0i)
              y7r = wn4r * (x2r - x2i)
              y7i = wn4r * (x2r + x2i)
              x0r = y1r + y5r
              x0i = y1i + y5i
              a(j + 2) = wk1r * x0r - wk1i * x0i
              a(j + 3) = wk1r * x0i + wk1i * x0r
              x0r = y1r - y5r
              x0i = y1i - y5i
              a(j + 10) = wk5r * x0r - wk5i * x0i
              a(j + 11) = wk5r * x0i + wk5i * x0r
              x0r = y3r - y7i
              x0i = y3i + y7r
              a(j + 6) = wk3r * x0r - wk3i * x0i
              a(j + 7) = wk3r * x0i + wk3i * x0r
              x0r = y3r + y7i
              x0i = y3i - y7r
              a(j + 14) = wk7r * x0r - wk7i * x0i
              a(j + 15) = wk7r * x0i + wk7i * x0r
              a(j) = y0r + y4r
              a(j + 1) = y0i + y4i
              x0r = y0r - y4r
              x0i = y0i - y4i
              a(j + 8) = wk4r * x0r - wk4i * x0i
              a(j + 9) = wk4r * x0i + wk4i * x0r
              x0r = y2r - y6i
              x0i = y2i + y6r
              a(j + 4) = wk2r * x0r - wk2i * x0i
              a(j + 5) = wk2r * x0i + wk2i * x0r
              x0r = y2r + y6i
              x0i = y2i - y6r
              a(j + 12) = wk6r * x0r - wk6i * x0i
              a(j + 13) = wk6r * x0i + wk6i * x0r
          end do
      end if
      end
!
      subroutine cftmdl(n, l, a, w)
      integer n, l, j, j1, j2, j3, j4, j5, j6, j7, k, k1, m
      real*8 a(0 : n - 1), w(0 : *)
      real*8 wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i
      real*8 wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      m = 8 * l
      wn4r = w(2)
      do j = 0, l - 2, 2
          j1 = j + l
          j2 = j1 + l
          j3 = j2 + l
          j4 = j3 + l
          j5 = j4 + l
          j6 = j5 + l
          j7 = j6 + l
          x0r = a(j) + a(j1)
          x0i = a(j + 1) + a(j1 + 1)
          x1r = a(j) - a(j1)
          x1i = a(j + 1) - a(j1 + 1)
          x2r = a(j2) + a(j3)
          x2i = a(j2 + 1) + a(j3 + 1)
          x3r = a(j2) - a(j3)
          x3i = a(j2 + 1) - a(j3 + 1)
          y0r = x0r + x2r
          y0i = x0i + x2i
          y2r = x0r - x2r
          y2i = x0i - x2i
          y1r = x1r - x3i
          y1i = x1i + x3r
          y3r = x1r + x3i
          y3i = x1i - x3r
          x0r = a(j4) + a(j5)
          x0i = a(j4 + 1) + a(j5 + 1)
          x1r = a(j4) - a(j5)
          x1i = a(j4 + 1) - a(j5 + 1)
          x2r = a(j6) + a(j7)
          x2i = a(j6 + 1) + a(j7 + 1)
          x3r = a(j6) - a(j7)
          x3i = a(j6 + 1) - a(j7 + 1)
          y4r = x0r + x2r
          y4i = x0i + x2i
          y6r = x0r - x2r
          y6i = x0i - x2i
          x0r = x1r - x3i
          x0i = x1i + x3r
          x2r = x1r + x3i
          x2i = x1i - x3r
          y5r = wn4r * (x0r - x0i)
          y5i = wn4r * (x0r + x0i)
          y7r = wn4r * (x2r - x2i)
          y7i = wn4r * (x2r + x2i)
          a(j1) = y1r + y5r
          a(j1 + 1) = y1i + y5i
          a(j5) = y1r - y5r
          a(j5 + 1) = y1i - y5i
          a(j3) = y3r - y7i
          a(j3 + 1) = y3i + y7r
          a(j7) = y3r + y7i
          a(j7 + 1) = y3i - y7r
          a(j) = y0r + y4r
          a(j + 1) = y0i + y4i
          a(j4) = y0r - y4r
          a(j4 + 1) = y0i - y4i
          a(j2) = y2r - y6i
          a(j2 + 1) = y2i + y6r
          a(j6) = y2r + y6i
          a(j6 + 1) = y2i - y6r
      end do
      if (m .lt. n) then
          wk1r = w(4)
          wk1i = w(5)
          do j = m, l + m - 2, 2
              j1 = j + l
              j2 = j1 + l
              j3 = j2 + l
              j4 = j3 + l
              j5 = j4 + l
              j6 = j5 + l
              j7 = j6 + l
              x0r = a(j) + a(j1)
              x0i = a(j + 1) + a(j1 + 1)
              x1r = a(j) - a(j1)
              x1i = a(j + 1) - a(j1 + 1)
              x2r = a(j2) + a(j3)
              x2i = a(j2 + 1) + a(j3 + 1)
              x3r = a(j2) - a(j3)
              x3i = a(j2 + 1) - a(j3 + 1)
              y0r = x0r + x2r
              y0i = x0i + x2i
              y2r = x0r - x2r
              y2i = x0i - x2i
              y1r = x1r - x3i
              y1i = x1i + x3r
              y3r = x1r + x3i
              y3i = x1i - x3r
              x0r = a(j4) + a(j5)
              x0i = a(j4 + 1) + a(j5 + 1)
              x1r = a(j4) - a(j5)
              x1i = a(j4 + 1) - a(j5 + 1)
              x2r = a(j6) + a(j7)
              x2i = a(j6 + 1) + a(j7 + 1)
              x3r = a(j6) - a(j7)
              x3i = a(j6 + 1) - a(j7 + 1)
              y4r = x0r + x2r
              y4i = x0i + x2i
              y6r = x0r - x2r
              y6i = x0i - x2i
              x0r = x1r - x3i
              x0i = x1i + x3r
              x2r = x1r + x3i
              x2i = x3r - x1i
              y5r = wk1i * x0r - wk1r * x0i
              y5i = wk1i * x0i + wk1r * x0r
              y7r = wk1r * x2r + wk1i * x2i
              y7i = wk1r * x2i - wk1i * x2r
              x0r = wk1r * y1r - wk1i * y1i
              x0i = wk1r * y1i + wk1i * y1r
              a(j1) = x0r + y5r
              a(j1 + 1) = x0i + y5i
              a(j5) = y5i - x0i
              a(j5 + 1) = x0r - y5r
              x0r = wk1i * y3r - wk1r * y3i
              x0i = wk1i * y3i + wk1r * y3r
              a(j3) = x0r - y7r
              a(j3 + 1) = x0i + y7i
              a(j7) = y7i - x0i
              a(j7 + 1) = x0r + y7r
              a(j) = y0r + y4r
              a(j + 1) = y0i + y4i
              a(j4) = y4i - y0i
              a(j4 + 1) = y0r - y4r
              x0r = y2r - y6i
              x0i = y2i + y6r
              a(j2) = wn4r * (x0r - x0i)
              a(j2 + 1) = wn4r * (x0i + x0r)
              x0r = y6r - y2i
              x0i = y2r + y6i
              a(j6) = wn4r * (x0r - x0i)
              a(j6 + 1) = wn4r * (x0i + x0r)
          end do
          k1 = 4
          do k = 2 * m, n - m, m
              k1 = k1 + 4
              wk1r = w(k1)
              wk1i = w(k1 + 1)
              wk2r = w(k1 + 2)
              wk2i = w(k1 + 3)
              wtmp = 2 * wk2i
              wk3r = wk1r - wtmp * wk1i
              wk3i = wtmp * wk1r - wk1i
              wk4r = 1 - wtmp * wk2i
              wk4i = wtmp * wk2r
              wtmp = 2 * wk4i
              wk5r = wk3r - wtmp * wk1i
              wk5i = wtmp * wk1r - wk3i
              wk6r = wk2r - wtmp * wk2i
              wk6i = wtmp * wk2r - wk2i
              wk7r = wk1r - wtmp * wk3i
              wk7i = wtmp * wk3r - wk1i
              do j = k, l + k - 2, 2
                  j1 = j + l
                  j2 = j1 + l
                  j3 = j2 + l
                  j4 = j3 + l
                  j5 = j4 + l
                  j6 = j5 + l
                  j7 = j6 + l
                  x0r = a(j) + a(j1)
                  x0i = a(j + 1) + a(j1 + 1)
                  x1r = a(j) - a(j1)
                  x1i = a(j + 1) - a(j1 + 1)
                  x2r = a(j2) + a(j3)
                  x2i = a(j2 + 1) + a(j3 + 1)
                  x3r = a(j2) - a(j3)
                  x3i = a(j2 + 1) - a(j3 + 1)
                  y0r = x0r + x2r
                  y0i = x0i + x2i
                  y2r = x0r - x2r
                  y2i = x0i - x2i
                  y1r = x1r - x3i
                  y1i = x1i + x3r
                  y3r = x1r + x3i
                  y3i = x1i - x3r
                  x0r = a(j4) + a(j5)
                  x0i = a(j4 + 1) + a(j5 + 1)
                  x1r = a(j4) - a(j5)
                  x1i = a(j4 + 1) - a(j5 + 1)
                  x2r = a(j6) + a(j7)
                  x2i = a(j6 + 1) + a(j7 + 1)
                  x3r = a(j6) - a(j7)
                  x3i = a(j6 + 1) - a(j7 + 1)
                  y4r = x0r + x2r
                  y4i = x0i + x2i
                  y6r = x0r - x2r
                  y6i = x0i - x2i
                  x0r = x1r - x3i
                  x0i = x1i + x3r
                  x2r = x1r + x3i
                  x2i = x1i - x3r
                  y5r = wn4r * (x0r - x0i)
                  y5i = wn4r * (x0r + x0i)
                  y7r = wn4r * (x2r - x2i)
                  y7i = wn4r * (x2r + x2i)
                  x0r = y1r + y5r
                  x0i = y1i + y5i
                  a(j1) = wk1r * x0r - wk1i * x0i
                  a(j1 + 1) = wk1r * x0i + wk1i * x0r
                  x0r = y1r - y5r
                  x0i = y1i - y5i
                  a(j5) = wk5r * x0r - wk5i * x0i
                  a(j5 + 1) = wk5r * x0i + wk5i * x0r
                  x0r = y3r - y7i
                  x0i = y3i + y7r
                  a(j3) = wk3r * x0r - wk3i * x0i
                  a(j3 + 1) = wk3r * x0i + wk3i * x0r
                  x0r = y3r + y7i
                  x0i = y3i - y7r
                  a(j7) = wk7r * x0r - wk7i * x0i
                  a(j7 + 1) = wk7r * x0i + wk7i * x0r
                  a(j) = y0r + y4r
                  a(j + 1) = y0i + y4i
                  x0r = y0r - y4r
                  x0i = y0i - y4i
                  a(j4) = wk4r * x0r - wk4i * x0i
                  a(j4 + 1) = wk4r * x0i + wk4i * x0r
                  x0r = y2r - y6i
                  x0i = y2i + y6r
                  a(j2) = wk2r * x0r - wk2i * x0i
                  a(j2 + 1) = wk2r * x0i + wk2i * x0r
                  x0r = y2r + y6i
                  x0i = y2i - y6r
                  a(j6) = wk6r * x0r - wk6i * x0i
                  a(j6 + 1) = wk6r * x0i + wk6i * x0r
              end do
          end do
      end if
      end
!
      subroutine rftfsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi
      m = n / 2
      ks = 2 * nc / m
      kk = 0
      do j = 2, m - 2, 2
          k = n - j
          kk = kk + ks
          wkr = 0.5d0 - c(nc - kk)
          wki = c(kk)
          xr = a(j) - a(k)
          xi = a(j + 1) + a(k + 1)
          yr = wkr * xr - wki * xi
          yi = wkr * xi + wki * xr
          a(j) = a(j) - yr
          a(j + 1) = a(j + 1) - yi
          a(k) = a(k) + yr
          a(k + 1) = a(k + 1) - yi
      end do
      end
!
      subroutine rftbsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi
      a(1) = -a(1)
      m = n / 2
      ks = 2 * nc / m
      kk = 0
      do j = 2, m - 2, 2
          k = n - j
          kk = kk + ks
          wkr = 0.5d0 - c(nc - kk)
          wki = c(kk)
          xr = a(j) - a(k)
          xi = a(j + 1) + a(k + 1)
          yr = wkr * xr + wki * xi
          yi = wkr * xi - wki * xr
          a(j) = a(j) - yr
          a(j + 1) = yi - a(j + 1)
          a(k) = a(k) + yr
          a(k + 1) = yi - a(k + 1)
      end do
      a(m + 1) = -a(m + 1)
      end
!
      subroutine dctsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
      m = n / 2
      ks = nc / n
      kk = 0
      do j = 1, m - 1
          k = n - j
          kk = kk + ks
          wkr = c(kk) - c(nc - kk)
          wki = c(kk) + c(nc - kk)
          xr = wki * a(j) - wkr * a(k)
          a(j) = wkr * a(j) + wki * a(k)
          a(k) = xr
      end do
      a(m) = c(0) * a(m)
      end
!
      subroutine dstsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
      m = n / 2
      ks = nc / n
      kk = 0
      do j = 1, m - 1
          k = n - j
          kk = kk + ks
          wkr = c(kk) - c(nc - kk)
          wki = c(kk) + c(nc - kk)
          xr = wki * a(k) - wkr * a(j)
          a(k) = wkr * a(k) + wki * a(j)
          a(j) = xr
      end do
      a(m) = c(0) * a(m)
      end
!

max_stars_repo_name: JunYangNXP/tensorflow

max_issues_repo_name: JunYangNXP/tensorflow

max_forks_repo_name: JunYangNXP/tensorflow


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE RMEW5JJ(J1,J2,K1,K2,COEF)
!                                                                  *
!   Written by  G. Gaigalas                                        *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE CONS_C,          ONLY: ZERO
      USE ribojj_C
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER,      INTENT(IN)  :: J1, J2, K1, K2
      REAL(DOUBLE), INTENT(OUT) :: COEF
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER                 :: JI1, JI2
      INTEGER, DIMENSION(6)   :: I00, I01
      INTEGER, DIMENSION(3,3) :: I12P,I12A,I03P, &
                       I03A,I14P,I14A,I05P,I05A
!-----------------------------------------------
      DATA I00/54,12,30,12,30,54/
      DATA I01/126,12,198,0,48,288/
      DATA I12P/420,360,270,360,2*0,-270,2*0/
      DATA I12A/0,1960,0,-1960,-1000,2430,0,2430,1980/
      DATA I03P/-882,3*0,384,-400,0,400,286/
      DATA I03A/4*0,162,-300,0,-300,22/
      DATA I14P/3780,-720,-4950,-720,2*0,4950,2*0/
      DATA I14A/2*0,3528,0,2430,1980,-3528,1980,-7722/
      DATA I05P/-1386,4*0,-440,0,440,-1430/
      DATA I05A/5*0,330,0,330,572/
!
      COEF=ZERO
      IF(IMPTJJ(J1) /= IMPTJJ(J2)) RETURN
      IF(J1 < 6 .OR. J1 > 11) RETURN
      IF(K1 == 0 .AND. K2 == 0) THEN
        IF(J1 /= J2) RETURN
        COEF=-DSQRT(DBLE(I00(J1-5)))
      ELSEIF(K1 == 1 .AND. K2 == 0) THEN
        IF(J1 /= J2) RETURN
        IF(J1 == 6) THEN
           COEF=-DSQRT(DBLE(48))
        ELSEIF(J1 == 9) THEN
           COEF=-DSQRT(DBLE(20))
        ELSEIF(J1 == 10) THEN
           COEF=-DSQRT(DBLE(10))
        ELSEIF(J1 == 11) THEN
           COEF=-DSQRT(DBLE(18))
        ENDIF
      ELSEIF(K1 == 0 .AND. K2 == 1) THEN
        IF(J1 /= J2) RETURN
        COEF=-DSQRT(DBLE(I01(J1-5))/DBLE(7))
      ELSEIF(K1 == 1 .AND. K2 == 2) THEN
        IF(J1 < 9) THEN
           JI1=J1-5
           JI2=J2-5
           IF(I12P(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I12P(JI1,JI2))/DBLE(7))
           ELSE
              COEF=-DSQRT(-DBLE(I12P(JI1,JI2))/DBLE(7))
           ENDIF
        ELSE
           JI1=J1-8
           JI2=J2-8
           IF(I12A(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I12A(JI1,JI2))/DBLE(49))
           ELSE
              COEF=-DSQRT(-DBLE(I12A(JI1,JI2))/DBLE(49))
           ENDIF
        ENDIF
      ELSEIF(K1 == 0 .AND. K2 == 3) THEN
        IF(J1 < 9) THEN
           JI1=J1-5
           JI2=J2-5
           IF(I03P(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I03P(JI1,JI2))/DBLE(21))
           ELSE
              COEF=-DSQRT(-DBLE(I03P(JI1,JI2))/DBLE(21))
           ENDIF
        ELSE
           JI1=J1-8
           JI2=J2-8
           IF(I03A(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I03A(JI1,JI2))/DBLE(7))
           ELSE
              COEF=-DSQRT(-DBLE(I03A(JI1,JI2))/DBLE(7))
           ENDIF
        ENDIF
      ELSEIF(K1 == 1 .AND. K2 == 4) THEN
        IF(J1 < 9) THEN
           JI1=J1-5
           JI2=J2-5
           IF(I14P(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I14P(JI1,JI2))/DBLE(35))
           ELSE
              COEF=-DSQRT(-DBLE(I14P(JI1,JI2))/DBLE(35))
           ENDIF
        ELSE
           JI1=J1-8
           JI2=J2-8
           IF(I14A(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I14A(JI1,JI2))/DBLE(49))
           ELSE
              COEF=-DSQRT(-DBLE(I14A(JI1,JI2))/DBLE(49))
           ENDIF
        ENDIF
      ELSEIF(K1 == 0 .AND. K2 == 5) THEN
        IF(J1 < 9) THEN
           JI1=J1-5
           JI2=J2-5
           IF(I05P(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I05P(JI1,JI2))/DBLE(21))
           ELSE
              COEF=-DSQRT(-DBLE(I05P(JI1,JI2))/DBLE(21))
           ENDIF
        ELSE
           JI1=J1-8
           JI2=J2-8
           IF(I05A(JI1,JI2) >= 0) THEN
              COEF=DSQRT(DBLE(I05A(JI1,JI2))/DBLE(7))
           ELSE
              COEF=-DSQRT(-DBLE(I05A(JI1,JI2))/DBLE(7))
           ENDIF
        ENDIF
      ELSE
         WRITE(0,'(A,4I5)') ' J1 J2 = ',J1,J2
         WRITE(0,'(A)') ' ERROR IN SUB. RMEW5JJ '
         STOP
      ENDIF
      RETURN
      END SUBROUTINE RMEW5JJ

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of Sandia Corporation nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: locout.f,v 1.1 1991/02/21 15:43:58 gdsjaar Exp $
C $Log: locout.f,v $
C Revision 1.1  1991/02/21 15:43:58  gdsjaar
C Initial revision
C
      SUBROUTINE LOCOUT (TYPE, NDIM, NODEL, TOLER, SORT, P1, P2, BOUND)
      DIMENSION P1(NDIM), P2(NDIM), TOLER(2)
      CHARACTER*(*) NODEL, BOUND, SORT, TYPE
      CHARACTER*16  BNAME
      include 'nu_io.blk'
C
      IF (NODEL(:1) .EQ. 'E') THEN
         NODEL = 'Elements'
      ELSE
         NODEL = 'Nodes'
      ENDIF
C
      IF (BOUND(:3) .EQ. 'BOU') THEN
         BNAME = 'Bounded Search'
      ELSE
         BNAME = 'Unbounded Search'
      END IF
C
      DO 10 IO=IOMIN, IOMAX
         WRITE (IO, 20) NODEL(:LENSTR(NODEL)), TOLER(1), TOLER(2),
     *      TYPE(:LENSTR(TYPE))
   20 FORMAT (//' Locating all ',A,' at a distance ',1PE15.8,
     *   ' plus/minus ',1PE15.8,/' from the ',A)
C
      IF (TYPE .EQ. 'LINE') THEN
         IF (NDIM .EQ. 2) THEN
            WRITE (IO, 30) (P1(I),I=1,NDIM), (P2(I),I=1,NDIM), BNAME
         ELSE
            WRITE (IO, 40) (P1(I),I=1,NDIM), (P2(I),I=1,NDIM), BNAME
         END IF
      ELSE IF (TYPE .EQ. 'POINT') THEN
         IF (NDIM .EQ. 2) THEN
            WRITE (IO, 50) (P1(I), I=1, NDIM)
         ELSE
            WRITE (IO, 60) (P1(I), I=1, NDIM)
         END IF
      ELSE IF (TYPE .EQ. 'PLANE') THEN
      A = P2(1)
      B = P2(2)
      C = P2(3)
      D = A * P1(1) + B * P1(2) + C * P1(3)
C
            WRITE (IO, 70) A, B, C, D
      END IF
      WRITE (IO, 80) SORT(:LENSTR(SORT))
   10 CONTINUE

   30 FORMAT (' from Point (',2(1PE11.3),')',/
     *        ' to   Point (',2(1PE11.3),')',2X,A)
   40 FORMAT (' from Point (',3(1PE11.3),')',/
     *        ' to   Point (',3(1PE11.3),')',2X,A)
   50 FORMAT (' (',2(1PE11.3),')')
   60 FORMAT (' (',3(1PE11.3),')')
   70 FORMAT (' ',1PE15.8,' X + ',1PE15.8,' Y + ',1PE15.8,
     *    ' Z = ',1PE15.8)
   80 FORMAT (' Sorted on field ',A,/)
C
      RETURN
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


! calculate bulk's energy band using wannier TB method
  subroutine ek_bulk2D

     use mpi
     use para

     implicit none

     integer :: ik, i, j, ib
     integer :: nkx, nky
	  integer :: knv3
     integer :: ierr
     integer :: nwann
     real(Dp) :: k(3)
     real(Dp) :: k11(3), k12(3)
     real(Dp) :: k21(3), k22(3)
     real(dp) :: sx
     real(dp) :: sy
     real(dp) :: sz
     real(Dp) :: W(Num_wann)
     real(dp) :: kxmin, kxmax, kymin, kymax
     real(dp), allocatable :: kxy(:,:)
     real(dp), allocatable :: kxy_shape(:,:)

     ! Hamiltonian of bulk system
     complex(Dp) :: Hamk_bulk(Num_wann,Num_wann)

     ! eigen value of H
	  real(dp), allocatable :: gap(:)
	  real(dp), allocatable :: gap_mpi(:)
	  real(dp), allocatable :: eigv(:,:)
	  real(dp), allocatable :: eigv_mpi(:,:)
	  real(dp), allocatable :: spin(:, :, :)
	  real(dp), allocatable :: spin_mpi(:, :, :)

     complex(dp), allocatable :: sigmax(:, :)
     complex(dp), allocatable :: sigmay(:, :)
     complex(dp), allocatable :: sigmaz(:, :)
     complex(dp), allocatable :: psi(:)

     nkx= Nk
     nky= Nk
     knv3= nkx*nky
     allocate( psi( Num_wann))
     allocate( sigmax( Num_wann, Num_wann))
     allocate( sigmay( Num_wann, Num_wann))
     allocate( sigmaz( Num_wann, Num_wann))
     allocate( spin    (3, Num_wann, knv3))
     allocate( spin_mpi(3, Num_wann, knv3))
     spin = 0d0
     spin_mpi = 0d0
     sigmax= 0d0
     sigmay= 0d0
     sigmaz= 0d0

     allocate( kxy(3, knv3))
     allocate( kxy_shape(3, knv3))
     allocate( gap     (   knv3))
     allocate( gap_mpi (   knv3))
     allocate( eigv    (4, knv3))
     allocate( eigv_mpi(4, knv3))
     kxy = 0d0
     gap    = 0d0
     gap_mpi= 0d0
	  eigv    = 0d0
	  eigv_mpi= 0d0

     nwann= Num_wann/2
     print *, 'nwann', nwann

     if (SOC==0) stop 'you should set soc=0 in the input file'
     !> construct spin matrix
     do i= 1, nwann
        sigmax(i, i+ nwann)= 1d0
        sigmax(i+ nwann, i)= 1d0
        sigmay(i, i+ nwann)= -zi
        sigmay(i+ nwann, i)= zi
        sigmaz(i, i)= 1d0
        sigmaz(i+ nwann, i+ nwann)= -1d0
     enddo

    !k11=(/-0.5d0,  0.0d0, -0.5d0/) !
    !k12=(/ 0.5d0,  0.0d0,  0.5d0/) ! X ! TB
     k11=(/-0.0d0,  0.0d0, -0.0d0/) !
     k12=(/ 1.0d0,  0.0d0,  1.0d0/) ! X ! TB
    !k12=(/ 0.0d0, -1.0d0, -1.0d0/) ! X ! DFT
     k21=(/-0.5d0, -0.5d0,  0.0d0/) !
     k22=(/ 0.5d0,  0.5d0,  0.0d0/) ! Z
    !k21=(/-0.0d0,  0.5d0,  0.5d0/) !
    !k22=(/ 0.0d0, -0.5d0, -0.5d0/) ! Y

     ik =0
     do i= 1, nkx
     do j= 1, nky
        ik =ik +1
        kxy(:, ik)= k11+(k12-k11)*(i-1)/dble(nkx-1)+  k21+ (k22-k21)*(j-1)/dble(nky-1)
        kxy_shape(:, ik)= kxy(1, ik)* Kua+ kxy(2, ik)* Kub+ kxy(3, ik)* Kuc
     enddo
     enddo

     do ik= 1+cpuid, knv3, num_cpu
	    !if (cpuid==0) print * , ik

        k = kxy(:, ik)

        ! calculation bulk hamiltonian
        Hamk_bulk= 0d0
        call ham_bulk(k, Hamk_bulk)

        !> diagonalization by call zheev in lapack
        W= 0d0
        call eigensystem_c( 'V', 'U', Num_wann ,Hamk_bulk, W)

        do ib= 1, Num_wann
           psi(:)= Hamk_bulk(:, ib)
          !print *, sum(conjg(psi(:))*psi(:))
           sx= 0d0
           sy= 0d0
           sz= 0d0
           do i= 1, Num_wann
              do j= 1, Num_wann
                 sx= sx+ conjg(psi(i))* sigmax(i, j)* psi(j)
                 sy= sy+ conjg(psi(i))* sigmay(i, j)* psi(j)
                 sz= sz+ conjg(psi(i))* sigmaz(i, j)* psi(j)
              enddo ! j
           enddo ! i
           spin(1, ib, ik)= sx
           spin(2, ib, ik)= sy
           spin(3, ib, ik)= sz
          !print *, W(ib),sx*sx+sy*sy+sz*sz
        enddo ! ib

        !> for 44 bands
        eigv(:, ik)= W(Numoccupied-1:Numoccupied+ 2)
        gap (   ik)= W(Numoccupied+1)- W(Numoccupied)
     enddo

     call mpi_allreduce(eigv,eigv_mpi,size(eigv),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(gap ,gap_mpi,size(gap ),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(spin ,spin_mpi,size(spin ),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)

     if (cpuid==0)then
        open(unit=14, file='bulkek2D.dat')
        do ik=1, knv3
           write(14, '(1000f19.9)')kxy_shape(:, ik), eigv_mpi(:, ik), (spin_mpi(:, ib, ik), ib=Numoccupied-1, Numoccupied+2)
        enddo
        close(14)
        open(unit=15, file='gap2D.dat')
        do ik=1, knv3
           write(15, '(1000f19.9)')kxy(:, ik), gap_mpi(ik)
        enddo
        close(15)
     endif

   return
   end subroutine ek_bulk2D

   !> using green's function
  subroutine ek_bulk2D_spin

     use mpi
     use para

     implicit none

     integer :: ik, i, j
     integer :: nkx, nky
	  integer :: nwann
	  integer :: knv3
     integer :: ierr
     real(Dp) :: k(3)
     real(dp) :: kx
     real(dp) :: ky
     real(dp) :: kz
     real(Dp) :: k11(3), k12(3)
     real(Dp) :: k21(3), k22(3)
     real(Dp) :: W(Num_wann)
     real(dp) :: kxmin, kxmax, kymin, kymax
     real(dp), allocatable :: kxy(:,:)
     real(dp), allocatable :: kxy_shape(:,:)

     ! Hamiltonian of bulk system
     complex(Dp) :: Hamk_bulk(Num_wann,Num_wann)

     ! eigen value of H
	  real(dp), allocatable :: gap(:)
	  real(dp), allocatable :: gap_mpi(:)
	  real(dp), allocatable :: eigv(:,:)
	  real(dp), allocatable :: eigv_mpi(:,:)
     real(dp), allocatable :: dos(:)
     real(dp), allocatable :: dos_mpi(:)
	  real(dp), allocatable :: spin(:, :)
	  real(dp), allocatable :: spin_mpi(:, :)

     complex(dp), allocatable :: sigmax(:, :)
     complex(dp), allocatable :: sigmay(:, :)
     complex(dp), allocatable :: sigmaz(:, :)
     complex(dp), allocatable :: psi(:)
     complex(dp), allocatable :: mat1(:, :)

     nkx= Nk
     nky= Nk
     knv3= nkx*nky
     allocate( psi( Num_wann))
     allocate( sigmax( Num_wann, Num_wann))
     allocate( sigmay( Num_wann, Num_wann))
     allocate( sigmaz( Num_wann, Num_wann))
     allocate( mat1  ( Num_wann, Num_wann))
     allocate( dos     (knv3))
     allocate( dos_mpi (knv3))
     allocate( spin    (3, knv3))
     allocate( spin_mpi(3, knv3))

     allocate( kxy(3, knv3))
     allocate( kxy_shape(3, knv3))
     allocate( gap     (   knv3))
     allocate( gap_mpi (   knv3))
     allocate( eigv    (4, knv3))
     allocate( eigv_mpi(4, knv3))
     kxy     = 0d0
     gap     = 0d0
     gap_mpi = 0d0
	  eigv    = 0d0
	  eigv_mpi= 0d0
     spin    = 0d0
     spin_mpi= 0d0
     dos    = 0d0
     dos_mpi= 0d0
     sigmax= 0d0
     sigmay= 0d0
     sigmaz= 0d0

     nwann= Num_wann/2

     if (SOC==0) stop 'you should set soc=0 in the input file'
     !> construct spin matrix
     do i= 1, nwann
        sigmax(i, i+ nwann)= 1d0
        sigmax(i+ nwann, i)= 1d0
        sigmay(i, i+ nwann)= -zi
        sigmay(i+ nwann, i)= zi
        sigmaz(i, i)= 1d0
        sigmaz(i+ nwann, i+ nwann)= -1d0
     enddo

     k11=(/-0.0d0,  0.0d0, -0.0d0/) !
     k12=(/ 1.0d0,  0.0d0,  1.0d0/) ! X
     k21=(/ 0.0d0,  0.0d0,  0.0d0/) !
     k22=(/ 1.0d0,  1.0d0,  0.0d0/) ! Z

     ik =0
     do i= 1, nkx
     do j= 1, nky
        ik =ik +1
        kxy(:, ik)= k11+(k12-k11)*(i-1)/dble(nkx-1)+  k21+ (k22-k21)*(j-1)/dble(nky-1)
        kxy_shape(:, ik)= kxy(1, ik)* Kua+ kxy(2, ik)* Kub+ kxy(3, ik)* Kuc
     enddo
     enddo

     do ik= 1+cpuid, nkx*nky, num_cpu
	     if (cpuid==0) print * , ik

        k = kxy(:, ik)

        ! calculation bulk hamiltonian
        Hamk_bulk= 0d0
        call ham_bulk(k, Hamk_bulk)

        !> diagonalization by call zheev in lapack
        Hamk_bulk= Hamk_bulk+ zi*eta_arc
        call inv(Num_wann, Hamk_bulk)
        do i=1, Num_wann
           dos(ik)= dos(ik) -aimag(Hamk_bulk(i, i))
        enddo

        call mat_mul(Num_wann, Hamk_bulk, sigmax, mat1)

        do i=1, Num_wann
           spin(1, ik)= spin(1, ik) -aimag(mat1(i, i))
        enddo

        call mat_mul(Num_wann, Hamk_bulk, sigmay, mat1)
        do i=1, Num_wann
           spin(2, ik)= spin(2, ik) -aimag(mat1(i, i))
        enddo

        call mat_mul(Num_wann, Hamk_bulk, sigmaz, mat1)
        do i=1, Num_wann
           spin(3, ik)= spin(3, ik) -aimag(mat1(i, i))
        enddo

        !> for 44 bands
        eigv(:, ik)= W(Numoccupied-1:Numoccupied+ 2)
        gap (   ik)= W(Numoccupied+1)- W(Numoccupied)
     enddo !ik

     call mpi_allreduce(eigv,eigv_mpi,size(eigv),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(gap ,gap_mpi,size(gap ),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(spin ,spin_mpi,size(spin ),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(dos ,dos_mpi,size(dos ),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)

     if (cpuid==0)then
        open(unit=14, file='bulkek2D.dat')
        do ik=1, knv3
           write(14, '(1000f20.6)')kxy_shape(:, ik), eigv_mpi(:, ik)
        enddo
        close(14)

        open(unit=15, file='gap2D.dat')
        do ik=1, knv3
           write(15, '(1000f20.6)')kxy_shape(:, ik), gap_mpi(ik)
        enddo
        close(15)

        open(unit=16, file='bulkspin2D.dat')
        ik= 0
        do i= 1, nkx
           do j= 1, nky
              ik= ik+ 1

              if (ik==1 .or. ik==knv3.or. dos_mpi(ik)<4d0) cycle
              if ((dos_mpi(ik)> dos_mpi(ik-1)) .and. (dos_mpi(ik)> dos_mpi(ik+1))) then
                !print *, ik, dos_mpi(ik)
                 write(16, '(10000f20.5)')kxy_shape(:, ik), dos_mpi(ik) , spin_mpi(:, ik)

              endif
           enddo
           write(16, '(a)')' '
        enddo
        close(16)

     endif

   return
   end subroutine ek_bulk2D_spin


text:
!    EmbryoMaker software (General Node Model)
!    Computational model to simulate morphogenetic processes in living organs and tissues.
!    Copyright (C) 2014 Miquel Marin-Riera, Miguel Brun-Usan, Roland Zimm, Tommi Välikangas & Isaac Salazar-Ciudad

!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    any later version.

!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.

!    You should have received a copy of the GNU General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.

module gnuplotter
implicit none
contains
subroutine  f2gp(nodoss,geness,imm)

integer :: n1,n2,j,ii,jj,im,ki,kii,rat
integer :: nodoss(10),geness(10)
integer :: plot_type                                 ! 1 for linear plot, 2 for log plot, 3 for log-log plot
 character(len=20) :: xlabel,ylabel                  ! plot axis labels and title
 character(len=40)  ::title1,title2,arxiv
!---------------

integer :: i,imm                                     ! number of lines i the plot points
integer :: ret

!---------------

xlabel='TIME'
ylabel='GENE CONCENTRATION'
plot_type=1                                         ! 1 for linear plot, 2 for log plot, 3 for log-log plot

!!!!!!!! it comes back if gnuplot is not installed
rat=0
call checkgnuplot(rat)
if(rat.eq.1)then
 write(*,*)'sorry, gene concentration can not be plotted because'
 write(*,*)'gnuplot is not installed'
 return
end if

! create gnuplot command file

OPEN(10,ACCESS='SEQUENTIAL',FILE='gp.txt')
write(10,*) 'set xlabel '//'"'//TRIM(xlabel)//'"'
write(10,*) 'set ylabel '//'"'//TRIM(ylabel)//'"'
write(10,*) 'set autoscale'
if (plot_type==2) write(10,*) 'set log y'
if (plot_type==3) then
   write(10,*) 'set log x' ;  write(10,*) 'set log y'
endif
write(10,*) 'set style line 1  lt 2 lw 1 pt 4 ps 2'

jj=0
do i=1,10
  ki=nodoss(i)
  if(ki.ne.0)then
    do j=1,10
      kii=geness(j)
      if(kii.ne.0)then
        jj=jj+1
        write(arxiv,*)"N",ki,"G",kii
        do im=1,40
          if (arxiv(im:im)==" ")then; arxiv(im:im)="_";end if
          arxiv(im:im)=arxiv(im:im)
        end do
        open(11111+jj,file=arxiv,status='unknown',action='read')
        write(title1,*)'node',ki,'gen',kii

        if(jj.lt.imm)then
          if(jj.eq.1)then
            write(10,*) 'plot "'//arxiv//'" using 1:2 with lines title "'//TRIM(title1)//'",\'
          else
            write(10,*) '"'//arxiv//'" using 1:2 with lines title "'//TRIM(title1)//'",\'
          end if
        else
          if(imm.gt.1)then
            write(10,*) '"'//arxiv//'" using 1:2 with lines title "'//TRIM(title1)//'"'
          else  ! only one gene in a node has been selected
            write(10,*) 'plot "'//arxiv//'" using 1:2 with lines title "'//TRIM(title1)//'"'
          end if
        end if
        close(11111+jj)
      end if
    end do
  end if
end do

 CLOSE(10,STATUS='KEEP')

! plot curve with gnuplot and cleanup files
ret=SYSTEM('gnuplot -persist gp.txt')
ret=SYSTEM('rm gp.txt')

end subroutine

!************************************************************************************************

subroutine checkgnuplot(rat)
 integer :: rat,ret,i
 character(len=10) :: output

 open(9999,file='check',status='unknown')
 ret=SYSTEM('dpkg -l | grep gnuplot > check')
 open(9959,file='check1',status='unknown')
 ret=SYSTEM('dpkg -l | du -bsh check > check1')
 read(9959,'(A)')output
 if(output(1:1).eq.'0')then
    rat=1
 end if
 ret=SYSTEM('rm check*')
! ret=SYSTEM('rm check1*')
end subroutine checkgnuplot

end module

max_issues_repo_name: sebjameswml/EmbryoMaker

max_forks_repo_name: sebjameswml/EmbryoMaker


text:
SUBROUTINE MULTIMODETRANSITIONAVG_C(D,NM,MODES_NUM,U_F_MODES,E_MULTIFLOQUET,D_BARE,U,INFO)
!!$   AVERAGE TIME EVOLUTION OPERATOR OF A MULTIMODE DRESSED SYSTEM. THE AVERAGE EVOLUTION OPERATOR IS WRITEN IN THE BASIS USED TO EXPRESS THE
!!$   MULTIMODE FLOQUET HAMILTONIAN
!!$   U : MATRIX OF AVERAGE TRANSITION PROBABILITIES
!!$
!!$  D              (IN)   : DIMENSION OF THE EXTENDED HILBERT SPACE (SIZE OF THE MULTIMODE FLOQUET MATRIX)
!!$  NM             (IN)   : NUMBER OF MODES
!!$  MODES_NUM      (IN)   : VECTOR (NM) INDICATING THE NUMBER OF HARMONICS OF EACH MODE
!!$  U_F_MODES      (IN)   : TRANSFORMATION, DIMENSOON (D,D)
!!$  E_MULTIFLOQUET (IN)   : MULTIMODE FLOQUET SPECTRUM
!!$  D_BARE         (IN)   : DIMENSION OF THE BARE HILBERT SPACE
!!$  U              (OUT)  :  MATRIX OF AVERAGE TRANSITION PROBABILITIES
!!$  INFO           (INOUT): (POSSIBLE) ERROR FLAG

  USE TYPES
  USE TYPES_C
  USE MODES_4F
  !  USE SUBINTERFACE_LAPACK

  IMPLICIT NONE
!  TYPE(MODE_C),DIMENSION(NM), INTENT(IN)     :: FIELD
  INTEGER,     DIMENSION(NM), INTENT(IN)     :: MODES_NUM

  INTEGER,                                    INTENT(IN)    :: D,D_BARE,NM ! DIMENSION OF THE MULTIMODE FLOQUET SPACE AND THE BARE BASIS
  INTEGER,                                    INTENT(INOUT) :: INFO
  DOUBLE PRECISION, DIMENSION(D),             INTENT(IN)    :: E_MULTIFLOQUET ! SET OF MULTIMODE FLOQUET ENERGIES, IN Hz, TO AVOID HBAR FACTORS
  COMPLEX*16,       DIMENSION(D,D),           INTENT(IN)    :: U_F_MODES   ! TRANFORMATION MATRIX BETWEEN DRESSED AND BARE BASIS
  DOUBLE PRECISION, DIMENSION(D_BARE,D_BARE), INTENT(OUT)   :: U           ! TIME AVERAGE TRANSITION MATRIX

  CALL MULTIMODETRANSITIONAVG(D,NM,COUPLING,MODES_NUM,U_F_MODES,E_MULTIFLOQUET,D_BARE,U,INFO)

END SUBROUTINE MULTIMODETRANSITIONAVG_C

!!$SUBROUTINE MULTIMODETRANSITIONPROBABILITY_DRESSEDBASIS(D_F,D,U_F,U_F_MODES,E_DRESSED,E_MULTIFLOQUET, &
!!$     & D_BARE,P,P2,INFO)
!!$
!!$  !D                : dimension of the multimode hilbert space
!!$  !D_F              : dimension of the dressed hilbert space
!!$  !U_F_MODES        : D x D matrix transformation to the multidressed eigenstates
!!$  !U_F              : D_F x D_F matrix transformation to the single-dressed eigenstates
!!$  !E_DRESSED        : single mode eigenvalues == dressed energies
!!$  !E_MULTIFLOQUET   : multimode eigenvalues
!!$  !D_BARE           : dimension of the bare hilbert space
!!$  !P                : matrix of Time and phase average transition probabilities
!!$  !INFO             : for info
!!$
!!$  ! The dressed basis is not always well defined.
!!$  !Time average transition probability
!!$
!!$  USE SUBINTERFACE_LAPACK
!!$!  USE FLOQUET
!!$  USE TYPES
!!$
!!$  IMPLICIT NONE
!!$  INTEGER,                                    INTENT(IN)     :: D_F,D,D_BARE
!!$  COMPLEX*16,       DIMENSION(D_F,D_F),       INTENT(INOUT)  :: U_F
!!$  COMPLEX*16,       DIMENSION(D,D),           INTENT(INOUT)  :: U_F_MODES
!!$  DOUBLE PRECISION, DIMENSION(D_F),           INTENT(IN)     :: E_DRESSED
!!$  DOUBLE PRECISION, DIMENSION(D),             INTENT(IN)     :: E_MULTIFLOQUET
!!$  !  COMPLEX*16,       DIMENSION(D_Bare,D_Bare),       INTENT(OUT)   :: P,P2
!!$  COMPLEX*16,       DIMENSION(D_F,D_F),       INTENT(OUT)    :: P,P2
!!$
!!$  INTEGER,                                    INTENT(INOUT)  :: INFO
!!$
!!$
!!$  TYPE(HARMONIC_FACTORS), DIMENSION(:),ALLOCATABLE :: U_n,U_MODES_n,U_nn,U_nn2
!!$  COMPLEX*16, DIMENSION(D,D)       :: U_AUX_F ! TEMPORAL VARIABLE, USED TO CHECK SUM_N U^n \DAGGER U^n = IDENTITY
!!$  !COMPLEX*16, DIMENSION(D_F,D_F)  :: U_AUX
!!$  COMPLEX*16, DIMENSION(D_F,D)     :: U_AUX
!!$  COMPLEX*16, DIMENSION(D,D_F)     :: U_AUX_2
!!$
!!$  DOUBLE PRECISION, DIMENSION(D)   :: PHASE_CONVENTION_MODES
!!$  DOUBLE PRECISION, DIMENSION(D_F) :: PHASE_CONVENTION_SINGLE
!!$
!!$  INTEGER n,m,index1,index2,i,j,index,index_,index3,index0
!!$  INTEGER DRESSED_HARMONICS, MULTIMODE_HARMONICS,D_AUX
!!$
!!$  INTEGER :: a,c,e,g,k,k_,n_,lambda,l,l_
!!$
!!$  DRESSED_HARMONICS   = D_F/D_BARE
!!$  MULTIMODE_HARMONICS =   D/D_BARE
!!$
!!$  ALLOCATE(U_n(DRESSED_HARMONICS))
!!$  ALLOCATE(U_MODES_n(MULTIMODE_HARMONICS))
!!$  ALLOCATE(U_nn(DRESSED_HARMONICS*MULTIMODE_HARMONICS))
!!$  ALLOCATE(U_nn2(DRESSED_HARMONICS*MULTIMODE_HARMONICS))
!!$
!!$  DO n=1,DRESSED_HARMONICS
!!$     ALLOCATE(U_n(n)%U(D_BARE,D_F))   ! All singlemode dressed states
!!$     !ALLOCATE(U_n(n)%U(D_BARE,D_BARE)) ! Select a manifold of singlemode dressed states
!!$     U_n(n)%U = 0.0
!!$  END DO
!!$
!!$  DO n=1,MULTIMODE_HARMONICS
!!$     ALLOCATE(U_MODES_n(n)%U(D_BARE,D)) !Full set of multimode dressed states
!!$     !ALLOCATE(U_MODES_n(n)%U(D_BARE,D_BARE)) !Single out a multimode dressed manifold
!!$     U_MODES_n(n)%U = 0.0
!!$     ALLOCATE(U_MODES_n(n)%n(MODES_NUM))
!!$     U_MODES_n(n)%n = 0
!!$  END DO
!!$
!!$  DO n=1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$     !ALLOCATE(U_nn(n)%U(D_F,D_BARE))
!!$     !U_nn(n)%U = 0.0
!!$     !ALLOCATE(U_nn2(n)%U(D_F,D_BARE))
!!$     !U_nn2(n)%U = 0.0
!!$     ALLOCATE(U_nn(n)%U(D_F,D))
!!$     U_nn(n)%U = 0.0
!!$     ALLOCATE(U_nn2(n)%U(D_F,D))
!!$     U_nn2(n)%U = 0.0
!!$  END DO
!!$  i  = 1
!!$  DO index1=1,DRESSED_HARMONICS
!!$     i = 1 + (index1-1)*D_BARE
!!$     U_n(index1)%U = U_F(i:i+D_BARE-1,1:D_F)
!!$     !index2 = D_BARE*N_FLOQUET(2)+1
!!$     !U_n(index1)%U = U_F(i:i+D_BARE-1,index2:index2+D_BARE-1)
!!$     !call write_matrix(ABS(U_n(index1)%U))
!!$     !
!!$  END DO
!!$
!!$
!!$  DO i=0,MULTIMODE_HARMONICS-1
!!$     U_MODES_n(i+1)%n = 0
!!$     index1 = i
!!$     DO j=2,MODES_NUM
!!$        U_MODES_n(i+1)%n(j)= -N_FLOQUET(j) + MOD(index1,(2*N_FLOQUET(j)+1))
!!$        index1 = INT(index1/(2*N_FLOQUET(j)+1))
!!$     END DO
!!$  END DO
!!$
!!$  i  = 1
!!$  DO index1=1,MULTIMODE_HARMONICS
!!$     i = 1 + (index1-1)*D_BARE
!!$     U_MODES_n(index1)%U = U_F_MODES(i:i+D_BARE-1,1:D)
!!$     !WRITE(*,*) INDEX1,U_MODES_n(index1)%n
!!$     !call write_matrix(ABS(U_MODES_n(index1)%U))
!!$
!!$     !index3 = MULTIMODE_HARMONICS*D_BARE/2.0 -1
!!$     !     !U_MODES_n(index1)%U = U_F_MODES(i:i+D_BARE-1,index3:index3+1)
!!$  END DO
!!$
!!$  index1 =1
!!$  DO n=1,DRESSED_HARMONICS
!!$     DO m=1,MULTIMODE_HARMONICS
!!$        U_nn(index1)%U  = ABS(MATMUL(TRANSPOSE(CONJG(U_n(n)%U)),U_MODES_n(m)%U))**2
!!$        !      CALL WRITE_MATRIX(ABS(U_n(n)%U))
!!$        !      CALL WRITE_MATRIX(ABS(U_MODES_n(m)%U))
!!$        !      CALL WRITE_MATRIX(ABS(U_nn(index1)%U))
!!$        index1 =  index1 + 1
!!$     END DO
!!$  END DO
!!$
!!$  !P = ABS(MATMUL(TRANSPOSE(CONJG(U_F)),U_F_MODES))**2
!!$  !P = MATMUL(P,TRANSPOSE(CONJG(P)))
!!$  !
!!$  INDEX0 = (DRESSED_HARMONICS -1)/2 + 1
!!$  INDEX1 = (MULTIMODE_HARMONICS -1)/2 + 1
!!$  !write(*,*) '#',index0,index1
!!$
!!$  U_AUX = ABS(MATMUL(TRANSPOSE(CONJG(U_n(INDEX0)%U)),U_MODES_n(INDEX1)%U))**2
!!$
!!$  !as in 10.10.17
!!$  P = 0.0
!!$  DO index1  = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS!,MULTIMODE_HARMONICS!
!!$     !write(*,*) size(U_aux,1),size(U_aux,2),size(u_nn(index1)%U,1),size(u_nn(index1)%U,2)
!!$     !P = P + MATMUL(U_AUX,transpose(conjg(U_nn(index1)%U)))
!!$     P = P + MATMUL(U_nn(index1)%U,TRANSPOSE(U_AUX))
!!$     !P = P + MATMUL(U_AUX,TRANSPOSE(U_AUX))
!!$  END DO
!!$
!!$
!!$  !  write(*,*) real(P(4,4)),real(P(4,5)),real(P(5,5))
!!$  !  P = 0.0
!!$  !  DO index1  = MULTIMODE_HARMONICS+1,MULTIMODE_HARMONICS*2!DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     P = P + MATMUL(U_nn(index1)%U,TRANSPOSE(U_AUX))
!!$  !  END DO
!!$  !  !CALL WRITE_MATRIX(ABS(P))
!!$  !
!!$  !  P = 0.0
!!$  !  DO index1  = 2*MULTIMODE_HARMONICS+1,MULTIMODE_HARMONICS*3!DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     P = P + MATMUL(U_nn(index1)%U,TRANSPOSE(U_AUX))
!!$  !  END DO
!!$  !  !CALL WRITE_MATRIX(ABS(P))
!!$  !
!!$  !
!!$  !  P = 0.0
!!$  !  DO index1  = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     P = P + MATMUL(U_nn(index1)%U,TRANSPOSE(U_AUX))
!!$  !  END DO
!!$  !  !CALL WRITE_MATRIX(ABS(P))
!!$  !
!!$  !  P = 0.0
!!$  !  P = MATMUL(U_AUX,TRANSPOSE(U_AUX))
!!$  !  !CALL WRITE_MATRIX(ABS(P))
!!$
!!$
!!$  ! as in 19.01.17 b.
!!$
!!$  !DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !   U_nn2(index1)%U = abs(U_nn(index1)%U)**2
!!$  !END DO
!!$  !P = 0.0
!!$  !DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !   DO index2 = 1,DRESSED_HARMONICS
!!$  !      INDEX3 = INDEX2 + (MULTIMODE_HARMONICS-1)/2
!!$  !      P =  P + MATMUL(U_nn2(index1)%U,TRANSPOSE(CONJG(U_nn2(index3)%U)))
!!$  !   END DO
!!$  !END DO
!!$
!!$
!!$  !! as in 17.01.17 ! selecting a dressed manifold
!!$
!!$  !  ALLOCATE(U_n(DRESSED_HARMONICS))
!!$  !  ALLOCATE(U_MODES_n(MULTIMODE_HARMONICS))
!!$  !  ALLOCATE(U_nn(DRESSED_HARMONICS*MULTIMODE_HARMONICS))
!!$  !  ALLOCATE(U_nn2(DRESSED_HARMONICS*MULTIMODE_HARMONICS))
!!$  !
!!$  !  DO n=1,DRESSED_HARMONICS
!!$  !     ALLOCATE(U_n(n)%U(D_BARE,D_BARE))
!!$  !     U_n(n)%U = 0.0
!!$  !  END DO
!!$  !
!!$  !  D_AUX = D_BARE
!!$  !  DO n=1,MULTIMODE_HARMONICS
!!$  !     ALLOCATE(U_MODES_n(n)%U(D_BARE,D_aux))
!!$  !     U_MODES_n(n)%U = 0.0
!!$  !     ALLOCATE(U_MODES_n(n)%n(MODES_NUM))
!!$  !     U_MODES_n(n)%n = 0
!!$  !  END DO
!!$  !
!!$  !  DO n=1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     ALLOCATE(U_nn(n)%U(D_BARE,D_AUX))
!!$  !     U_nn(n)%U = 0.0
!!$  !     ALLOCATE(U_nn2(n)%U(D_BARE,D_AUX))
!!$  !     U_nn2(n)%U = 0.0
!!$  !  END DO
!!$  !
!!$  !  DO i=0,MULTIMODE_HARMONICS-1
!!$  !     U_MODES_n(i+1)%n = 0
!!$  !     index1 = i
!!$  !     DO j=2,MODES_NUM
!!$  !        U_MODES_n(i+1)%n(j)= -N_FLOQUET(j) + MOD(index1,(2*N_FLOQUET(j)+1))
!!$  !        index1 = INT(index1/(2*N_FLOQUET(j)+1))
!!$  !     END DO
!!$  !  END DO
!!$  !
!!$  !
!!$  !  i  = 1
!!$  !  DO index1=1,DRESSED_HARMONICS
!!$  !     i =1 + (index1-1)*D_BARE
!!$  !     index2 = D_BARE*N_FLOQUET(2)+2
!!$  !     U_n(index1)%U = U_F(i:i+D_BARE-1,index2:index2+D_BARE-1)
!!$  !  END DO
!!$  !  i  = 1
!!$  !  DO index1=1,MULTIMODE_HARMONICS
!!$  !     i = 1 + (index1-1)*D_BARE
!!$  !     !U_MODES_n(index1)%U = U_F_MODES(i:i+D_BARE-1,1:D)
!!$  !     index = MULTIMODE_HARMONICS*D_BARE/2.0 + D_BARE+1
!!$  !     U_MODES_n(index1)%U = U_F_MODES(i:i+D_BARE-1,index:index+D_AUX-1)
!!$  !     if(index1.eq.5) then
!!$  !         P = U_MODES_n(index1)%U
!!$  !     end if
!!$  !  END DO
!!$  !
!!$  !  index1 =1
!!$  !  DO n=1,DRESSED_HARMONICS
!!$  !     DO m=1,MULTIMODE_HARMONICS
!!$  !        U_nn(index1)%U  = MATMUL(TRANSPOSE(CONJG(U_n(n)%U)),U_MODES_n(m)%U)
!!$  !        !write(*,*) n,m,index1
!!$  !        index1 =  index1 + 1
!!$  !     END DO
!!$  !  END DO
!!$  !
!!$  !!17.01.17  expressions where the sums are over a single set of multimode dressed states, page 1
!!$  !
!!$  !  DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     U_nn2(index1)%U = ABS(U_nn(index1)%U)**2
!!$  !!    CALL WRITE_MATRIX(ABS(U_nn2(index1)%U))
!!$  !  END DO
!!$  !
!!$  !  P = 0.0
!!$  !  DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     DO index2 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !              P      = P+MATMUL(U_nn2(INDEX1)%U,TRANSPOSE(U_nn2(INDEX2)%U))
!!$  ! !       CALL WRITE_MATRIX(REAL(P))
!!$  !     END DO
!!$  !  END DO
!!$
!!$  !13.01.17  expressions where the sums are over a single set of multimode dressed states, page 12
!!$
!!$  !  DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !     U_nn2(index1)%U = MATMUL(U_nn(index1)%U,TRANSPOSE(CONJG(U_nn(index1)%U)))
!!$  !  END DO
!!$  !
!!$  !  P = 0.0
!!$  !  aDO i=1,D_BARE
!!$  !     DO j=1,D_BARE
!!$  !        DO index1 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !           DO index2 = 1,DRESSED_HARMONICS*MULTIMODE_HARMONICS
!!$  !              P(i,j) = P(i,j) + U_nn2(index1)%U(j,j)*U_nn2(index2)%U(i,i)
!!$  !           END DO
!!$  !        END DO
!!$  !     END DO
!!$  !  END DO
!!$  !
!!$
!!$
!!$END SUBROUTINE MULTIMODETRANSITIONPROBABILITY_DRESSEDBASIS
!!$

max_stars_repo_name: gsinuco/OPENMMF

max_issues_repo_name: gsinuco/OPENMMF

max_forks_repo_name: gsinuco/OPENMMF


text:
! @@name:	tasking.11f
! @@type:	F-free
! @@compilable:	yes
! @@linkable:	no
! @@expect:	success
subroutine foo()
  integer :: x
  x = 2
!$omp task shared(x) mergeable
  x = x + 1
!$omp end task
!$omp taskwait
  print *, x     ! prints 3
end subroutine

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
C
C
C
C     ---------------------------------------------------
      SUBROUTINE GETTRIGGER(IERRT)
C     ---------------------------------------------------
      implicit none
      integer GFLT_var(20)
C     ---------------------------------------------------
C
C  COMMENTS FROM MIKE Wodarczyk
C
C  get_flt_data fills the array GFLT_VAR with the variables necessary
C    to apply restrictions to GFLT slots 46 and 47
C    As a side effect, the common block o1cpdt is filled with ALL GFLT
C    trigger quantities comming from the data (or monte carlo zgana).
C
C  Output:
C     integer GFLT_VAR(20)  ! array of GFLT quantities needed
C               GFLT_VAR(1) = o1evnt_Subtrg(1) bits 0 to 31   (Needed)
C               GFLT_VAR(2) = o1evnt_Subtrg(2) bite 32 to 63  (Needed)
C               GFLT_VAR(3) = o1subt_Subtrg(1) bits 0 to 31   (Not Needed)
C               GFLT_VAR(4) = o1subt_Subtrg(2) bite 32 to 63  (Not Needed)
C               GFLT_VAR(5) = AH_SRTD_s0_ti  srtd timing 0    (Not Needed)
C               GFLT_VAR(6) = AH_SRTD_s1_ti  srtd timing 1    (Not Needed)
C               GFLT_VAR(7) = AH_SRTD_s2_ti  srtd timing 2    (Not Needed)
C               GFLT_VAR(8) = AH_SRTD_s3_ti  srtd timing 3    (Not Needed)
C               GFLT_VAR(9) = RCAL_EMC_E     CFLT value        (Needed)
C               GFLT_VAR(10) = CAL_E         CFLT E total      (Needed)
C               GFLT_VAR(11) = Rcal_EMCTh    Rcal EMC Threshold bit energy   (Needed)
C               GFLT_VAR(12) = BCAL_EMC_E    Bcal Emc energy   (Needed)
C               GFLT_VAR(13) = CAL_Et        CFLT Et           (Needed)
C               GFLT_VAR(14) = RCAL_ISOE     Sum of Rcal isoe's   (Needed)
C               GFLT_VAR(15:20) = [ not used yet ]
C
C  Implicit output:
C      If you want other GFLT quantities for personal use then
C      place the following includes in your zuanal routine to
C      access other GFLT variables.
C
C  #include "o1cmn/gflcpd.inc"
C  C    integer o1cpdt(1000) is the array of ALL GFLT Variables
C  #include "FLT/GFLT/gfltpr.inc"
C  C    parameters in gfltpr allow you to index the array o1cpdt as follows
C  C    REMC = o1cpdt(RCAL_EMC_E)
C
C ======================================================================

#include "common.inc"
c#include "comnt.inc"
#include "fmckin.inc"
#include "partap.inc"
!  /*       <-- this is the necessary ADAMO definitions  */
#include "zrevt.inc"
!  /*        <-- event information  */
#include "zdskey.inc"
#include "gabffw.inc"
#include "vctrhl.inc"
#include "lmresu.inc"
#include "lmposd.inc"
#include "lmeb.inc"
#include "gfltpr.inc"
#include "gfltcr.inc"
#include "ccgsum.inc"
C#include "datatyp.inc"
#include "o1evnt.inc"
!  /*               GFLT Trigger information  */
#include "o1subt.inc"
!  /*               GFLT Subtriggers that fired  */
#include "o2sdec.inc"
#include "o2sbor.inc"
#include "o4sbor.inc"
#include "tltevt.inc"
C#include "gflpar.inc"
c#include "FLT/GFLT/gfltpr.inc"
!  /*               indices for GFLT common block words */
c#include "o1cmn/gflcpd.inc"
#include "gflcpd.inc"
!  /*               for O1CPDT  from o1recon */
c#include "FLT/GFLT/gfltcr.inc"
!  /*               for O1RCEX  */

c Values used by Adi but not be me (yet)
c      integer FLT(4), SLT, TLT(15), DST
      integer FLT(4), SLT, DST
       integer Ndata,GFLT_Array(999),bcn
c       integer RCAL_EMC_EBP,RCAL_EMC_E,
c     +         RCAL07_ISOL_E
c       PARAMETER ( RCAL07_ISOL_E = 101 )
c       PARAMETER ( RCAL_EMC_E = 96 )
c       PARAMETER ( RCAL_EMC_EBP = 95 )
c       INTEGER RCAL_EBP_STB
c            PARAMETER ( RCAL_EBP_STB = 92 )
c       INTEGER RCAL_E_OUT_BP
c            PARAMETER ( RCAL_E_OUT_BP = 93 )
c      INTEGER RCAL_E_AWAY_B
c            PARAMETER ( RCAL_E_AWAY_B = 94 )
c       INTEGER LUMI_EGAMMA
c            PARAMETER ( LUMI_EGAMMA = 159 )
c       INTEGER LUMI_EE
c            PARAMETER ( LUMI_EE = 160 )

      logical trigmontecarlo
      integer it,ierrt

      LOGICAL FIRST
      data first /.true./

      ierrt=0

      IF (first) then
       WRITE (*,*) ' ************* '
       WRITE (*,*) '   GETTRIGGER  '
       WRITE (*,*) ' ************* '
       first = .false.
      ENDIF

C>>>>>>>
C
C   Check if data or monte carlo
C
      if (coutab(FmcKin).gt.0) then
         trigmontecarlo = .TRUE.
      else
         trigmontecarlo = .FALSE.
      endif

C Trigger stuff
C
C
C First level trigger
      call vzero(flt,4)
      if(coutab(O1EVNT).gt.0)then
         call FETTAB(O1EVNT,ID,1)     !!! prescaled
         FLT(1) = O1EVNT_Subtrg(1)
         FLT(2) = O1EVNT_Subtrg(2)
      endif
      if(coutab(O1SUBT).gt.0)then
         call FETTAB(O1SUBT,ID,1)     !!! without prescales
         FLT(3) = O1SUBT_Subtrg(1)
         FLT(4) = O1SUBT_Subtrg(2)
      endif
C
C get CFLT data
C
C      call getfltd(imc, cale, bcal_emc, rcal_emc, remcth,
C     &     et, emc_e, tlumi_e, tlumi_g)
C
C      write(6, '(A, I3, 8I8)') 'data',
C     &     imc, cale, bcal_emc, rcal_emc, remcth, et, emc_e,
C     &     tlumi_e, tlumi_g
C
C
C Second level trigger
      SLT = 0
      IF ( COUTAB(O2SDEC).GT.0 ) THEN
         DO IT=1,COUTAB(O2SDEC)
            CALL FETTAB(O2SDEC,ID,IT)
            IF ( O2SDEC_SUBTRIGNO.EQ.2 )
     &           SLT = SLT + IBITS( O2SDEC_SUBTRIGTYPE(1),0,8 )
            IF ( O2SDEC_SUBTRIGNO.EQ.3 ) THEN
               IF ( BTEST(O2SDEC_SUBTRIGTYPE(1),5) ) SLT = SLT + 256
               IF ( BTEST(O2SDEC_SUBTRIGTYPE(1),6) ) SLT = SLT + 512
               IF ( BTEST(O2SDEC_SUBTRIGTYPE(1),7) ) SLT = SLT + 1024
               IF ( BTEST(O2SDEC_SUBTRIGTYPE(1),8) ) SLT = SLT + 2048
            END IF
         END DO
      END IF
C
C
C
C... Third level trigger
      call vzero(TLT,15)
      if(coutab(TLTEVT).gt.0)then
         call FETTAB(TLTEVT,ID,1)
C         TLT(4) = TLTEVT_Subtrg(4)
C         TLT(1) = TLTEVT_Subtrg(1)
         do iT=1,15
          TLT(IT) = TLTEVT_Subtrg(IT)
         enddo
      endif
C
C... DST bits
      dst = zdskey_tstam11
C
C
C
C.. additional trigger info
C -----------------------------------------------------------------
C  First get the GFLT values for things like RCAL_EMC_E and CAL_E.
C
         If (.not.trigMonteCarlo) then
            CALL O1GTCV(ierrt, O1CPDT)
         else
            call o1vzst(-1)
            call o1ezin(ierrt)
            if (ierrt .ne. 0) then
               write(6, *) 'ERROR getfltd: MC FLT data is missing'
            endif
C     First, just the dead-copy
            MAX_O1CPDT = MAX_O1RCEX
            TOT_O1CPDT = TOT_O1RCEX
            DO  IT = 1, MAX_O1RCEX
               O1CPDT(IT) = O1RCEX(IT)
            ENDDO
            DO  IT = MAX_O1RCEX+1, LIMIT_O1CPDT
               O1CPDT(IT) = -1000
            ENDDO
         endif

C
C
C
C
C  Now the common block O1CPDT is filled with GFLT information
C     the indices are given in the file "FLT/GFLT/gfltpr.inc"
C
C
C  There are tons of variables that you might want to save, but
C    I wrote this specifically to collect the variables used
C    for GFLT(46) in 1996 and 1997.  I will pack them in an
C    array suitable for passing to the routines that will
C    restrict GFLT(46) to only the parts used for a High Q2
C    F2 analysis.
C
C

         GFLT_VAR(5) = o1cpdt(AH_SRTD_s0_ti)
         GFLT_VAR(6) = o1cpdt(AH_SRTD_s1_ti)
         GFLT_VAR(7) = o1cpdt(AH_SRTD_s2_ti)
         GFLT_VAR(8) = o1cpdt(AH_SRTD_s3_ti)

         GFLT_VAR(9) = o1cpdt(RCAL_EMC_E)
         GFLT_VAR(10) = o1cpdt(CAL_E)
         GFLT_VAR(11) = o1cpdt(RCAL_EBP_STB) + o1cpdt(RCAL_E_OUT_BP)
         GFLT_VAR(12) = o1cpdt(BCAL_EMC_E)
         GFLT_VAR(13) = o1cpdt(CAL_ET)
         GFLT_VAR(14) = o1cpdt(RCAL04_ISOL_E)+o1cpdt(RCAL05_ISOL_E)+
     &        o1cpdt(RCAL06_ISOL_E)+o1cpdt(RCAL07_ISOL_E)

C -----------------------------------------------------------------
C
C
C
C
      call VZERO(TrigDat,15)
      call O1GTCV(Ndata,GFLT_Array)
      TrigDat(1) = Float(GFLT_Array(LUMI_EGAMMA))
      TrigDat(2) = Float(GFLT_Array(LUMI_EE))
C
      Trigdat(3) = Float(GFLT_Array(Rcal_EMC_E))
C
C      TrigDat(2) = Float(GFLT_Array(Rcal_EMC_EBP))
C      TrigDat(3) = Float(GFLT_Array(Rcal_E_OUT_BP))
C      TrigDat(4) = Float(GFLT_Array(Rcal_EBP_STB))
C
      if(coutab(CCGSUM).gt.0) then
         call FETTAB(CCGSUM,ID,1)
         TrigDat(4) = (CCGSUM_ETotalEmc + CCGSUM_ETotalHac)
     +              - (CCGSUM_PzEmc + CCGSUM_PzHac)
      endif
      if(coutab(TLTEVT).gt.0)then
         call FETTAB(TLTEVT,ID,1)
         TrigDat(5) = TLTEVT_Eminpz
      endif
      TrigDat(6) = Float(GFLT_VAR(5))
      TrigDat(7) = Float(GFLT_VAR(6))
      TrigDat(8) = Float(GFLT_VAR(7))
      TrigDat(9) = Float(GFLT_VAR(8))
      TrigDat(10) = Float(GFLT_VAR(9))
      TrigDat(11) = Float(GFLT_VAR(10))
      TrigDat(12) = Float(GFLT_VAR(11))
      TrigDat(13) = Float(GFLT_VAR(12))
      TrigDat(14) = Float(GFLT_VAR(13))
      TrigDat(15) = Float(GFLT_VAR(14))

C>>>>>>>

      RETURN
      END
C

ext: fpp


text:
!** Copyright (c) 2000, NVIDIA CORPORATION.  All rights reserved.
!**
!** Licensed under the Apache License, Version 2.0 (the "License");
!** you may not use this file except in compliance with the License.
!** You may obtain a copy of the License at
!**
!**     http://www.apache.org/licenses/LICENSE-2.0
!**
!** Unless required by applicable law or agreed to in writing, software
!** distributed under the License is distributed on an "AS IS" BASIS,
!** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!** See the License for the specific language governing permissions and
!** limitations under the License.
!
!       Modify test bs04.f90 to use NEWUNIT

	program p
	 implicit none
	 character*50 a(10),b(10)
	 integer ::i,errors=0,expect=0
	 integer*1 i1(10),j1
	 integer*2 i2(10),j2
	 integer*4 i4(10),j4
	 integer*8 i8(10),j8
	 integer*1 u1(4)
	 integer*2 u2(4)
	 integer*4 u3(4)
	 integer*8 u4(4)

	 a(1) = 'abc abc abc abc abc'
	 a(2) = 'abc def abc abc abc'
	 a(3) = 'abc def ghi abc abc'
	 a(4) = 'abc def ghi jkl abc'
	 a(5) = 'abc def ghi jkl mno'
	 a(6) = 'abc def ghi        '
	 a(7) = 'abc def            '
	 a(8) = 'abc                '
	 a(9) = '0000               '
	 a(10)= '-------------------'

	 open(newunit=u1(2),file='test.txt',access='DIRECT',&
		form='UNFORMATTED',recl=50)
	 do i = 1,10
	  write(u1(2),rec=i) a(i)
	 enddo
	 close(unit=u1(2))

	 i1 = 66
	 i1(8) = 50
	 i1(6) = 11
	 open(newunit=u2(2),file='test.txt',access='DIRECT',&
		action='READ',form='UNFORMATTED',recl=i1(8))
	 do j1 = 1,10
	  i1(3) = j1
	  read(u2(2),rec=i1(3)) b(i1(3))
	  if( b(i1(3)) .ne. a(i1(3)) ) then
	   print *,'i1 error ',i1(3),' is ',b(i1(3))
	   errors = errors + 1
	  endif
	 enddo
	 close(unit=u2(2))

	 i2 = 6666
	 i2(4) = 50
	 i2(9) = 12
	 open(newunit=u3(3),file='test.txt',access='DIRECT',&
		action='READ',form='UNFORMATTED',recl=i2(4))
	 do j2 = 1,10
	  i2(2) = j2
	  read(u3(3),rec=i2(2)) b(i2(2))
	  if( b(i2(2)) .ne. a(i2(2)) ) then
	   print *,'i2 error ',i2(2),' is ',b(i2(2))
	   errors = errors + 1
	  endif
	 enddo
	 close(unit=u3(3))

	 i4 = 6666666
	 i4(1) = 50
	 i4(3) = 14
	 open(newunit=u4(4),file='test.txt',access='DIRECT',&
		action='READ',form='UNFORMATTED',recl=i4(1))
	 do j4 = 1,10
	  i4(5) = j4
	  read(u4(4),rec=i4(5)) b(i4(5))
	  if( b(i4(5)) .ne. a(i4(5)) ) then
	   print *,'i4 error ',i4(5),' is ',b(i4(5))
	   errors = errors + 1
	  endif
	 enddo
	 close(unit=u4(4))

	 i8 = 6666666666
	 i8(9) = 50
	 i8(8) = 18
	 open(newunit=u1(1),file='test.txt',access='DIRECT',&
		action='READ',form='UNFORMATTED',recl=i8(9))
	 do j8 = 1,10
	  i8(7) = j8
	  read(u1(1),rec=i8(7)) b(i8(7))
	  if( b(i8(7)) .ne. a(i8(7)) ) then
	   print *,'i8 error ',i8(7),' is ',b(i8(7))
	   errors = errors + 1
	  endif
	 enddo
	 close(unit=u1(1))

	 !print *,errors,' errors'
	 call check(errors,expect,1)
	end


text:
subroutine units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2)
  use amr_commons
  use hydro_commons
  use cooling_module
  real(dp)::scale_nH,scale_T2,scale_t,scale_v,scale_d,scale_l
  !-----------------------------------------------------------------------
  ! Conversion factors from user units into cgs units
  ! For gravity runs, make sure that G=1 in user units.
  !-----------------------------------------------------------------------

  ! scale_d converts mass density from user units into g/cc
  scale_d = 3.9234d-18

  ! scale_t converts time from user units into seconds
  scale_t = 1.0/sqrt(6.67d-8*scale_d)

  ! scale_l converts distance from user units into cm
  scale_l = 4.947d+16

  ! scale_v convert velocity in user units into cm/s
  scale_v = scale_l / scale_t

  ! scale_T2 converts (P/rho) in user unit into (T/mu) in Kelvin
  scale_T2 = mH/kB * scale_v**2

  ! scale_nH converts rho in user units into nH in H/cc
  scale_nH = X/mH * scale_d

end subroutine units

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
Describe Users/Kicker here.

20080630 20:01:14 nbsp  Welcome to the Wiki. Please consider the importance of using your RealName. Your comments on Davis Swim and Fitness Club caused me to take a look at your http://daviswiki.org/Users/Kicker?actionuserinfo edit history. Users/JasonAller

20080712 23:25:18 nbsp Hi Kicker, I think your name is Scott. Could be Casey, or who ever.... you should use your real name. I hope you are happy with the comments you have made about our facility. I am not sure where all this anger is coming from but it would be nice to talk to you face to face instead of hidding behind a computer and trying to damage my business. You know I am a good person and depend on this business to take care of my family, do you have any children? If you have anything to say, I would appreciate you being honest with yourself and just look at what will become of spreading negative (untrue) comments about me personal comments that are hurtful and very ugly. I would appreciate if you would go to Peak Proformance and stop with the lashing, we get it...  Users/TammyHengel

20080713 10:09:59 nbsp Please read Welcome to the Wiki/Business Owner and then take a look at your behavior and see if your persistence about Davis Swim and Fitness Club is helping out, or causing problems. With free speech comes responsibility; please be responsible about your speech here. Other users have been BannedGroup/Talk banned because they could not be responsible about their behavior and speech. Users/JasonAller

20091116 17:47:47 nbsp November 16, 2009

Davis Swim and Fit
Attn: Tammy Hengel
Club Manager
Davis CA

Per your request here is documentation in writing that I am not now nor have I ever been a member of Davis Swim and Fit. I signed up to be a member of 24 Hr Fitness and never used or visited your facility. I am told that 24 Hr fitness sold my membership to your club in October of 2006 but continued to honor my membership at their clubs through early 2008.

Davis Swim and Fit was taking electronic fund  transfers from my Golden 1 Checking Account through September of 2009.  I was told by 24hr club managers that Davis Swim and Fit was a 24 hr Fitness Club that I was a 24hr Fitness member and that the debits were going to the club of original record.

I was never a member of Davis Swim and Fit I signed up to be a member of 24 hr Fitness and have never visited your club or even know where it is located.

I have turned this matter over to the Fraud dept at Golden 1 and will make every attempt to recover the $480.00 that was illegally taken from my account by your club.
Please consider this notice that I do not now nor ever want to be a member of Davis Swim and Fit and further attempts to withdraw funds from my account will not be processed. I demand a full refund for funds taken from October 2007 September 2008.

Most Sincerely,

Tamara Morgan
 Users/Sacramento24

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search