Class at_module In: src/at_module.f90

Included Modules

dc_message lumatrix

Public Instance methods

 a_Int_ag : real(8), dimension(size(ag,1)) ag : real(8), dimension(:,0:), intent(in)

————— 積分計算 ——————

[Source]

```    function a_Int_ag(ag)
real(8), dimension(:,0:), intent(in)     :: ag
real(8), dimension(size(ag,1))           :: a_Int_ag
integer :: i

if ( size(ag,2) < im+1 ) then
call MessageNotify('E','ae_ag',               'The Grid points of input data too small.')
elseif ( size(ag,2) > im+1 ) then
call MessageNotify('W','ae_ag',               'The Grid points of input data too large.')
endif

a_Int_ag = 0.0d0
do i=0,im
a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i)
enddo
end function a_Int_ag
```
 ag_at : double precision, dimension(size(at_data,1),0:im) at_data : double precision, dimension(:,:), intent(in)

—- 逆変換 —-

[Source]

```  function ag_at(at_data)  ! スペクトル -> 台形格子(2 次元)

double precision, dimension(:,:), intent(in)      :: at_data
double precision, dimension(size(at_data,1),0:im)  :: ag_at

real(8), dimension(size(at_data,1)*im) :: y
integer :: m

m = size(at_data,1)
if ( size(at_data,2)-1 < km ) then
call MessageNotify('E','ag_at',             'The Chebyshev dimension of input data too small.')
elseif ( size(at_data,2)-1 > km ) then
call MessageNotify('W','ag_at',             'The Chebyshev dimension of input data too large.')
endif

ag_at = 0.0D0
ag_at(:,0:km)=at_data
call fttctb(m,im,ag_at,y,it,t)

end function ag_at
```
 at_Dx_at : real(8), dimension(size(at_data,1),0:size(at_data,2)-1) at_data : real(8), dimension(:,0:), intent(in)

—- 微分計算 —-

[Source]

```  function at_Dx_at(at_data)
real(8), dimension(:,0:), intent(in)                    :: at_data
real(8), dimension(size(at_data,1),0:size(at_data,2)-1) :: at_Dx_at

integer :: m, k
integer :: nm, kmax

nm=size(at_data,1)
kmax=size(at_data,2)-1
if ( kmax  < km ) then
call MessageNotify('W','at_Dx_at',             'The Chebyshev dimension of input data too small.')
elseif ( kmax > km ) then
call MessageNotify('E','at_Dx_at',             'The Chebyshev dimension of input data too large.')
endif

if ( kmax == im ) then
do m=1,nm
at_Dx_at(m,kmax)   = 0.
at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax) /2
enddo
else
do m=1,nm
at_Dx_at(m,kmax)   = 0.
at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax)
! スタートはグリッド対応最大波数未満. Factor 1/2 不要
enddo
endif

do k=kmax-2,0,-1
do m=1,nm
at_Dx_at(m,k) = at_Dx_at(m,k+2) + 2*(k+1)*at_data(m,k+1)
enddo
enddo

do k=0,kmax
do m=1,nm
at_Dx_at(m,k) = 2/xl * at_Dx_at(m,k)
enddo
enddo

end function at_Dx_at
```
 i : integer, intent(in) k : integer, intent(in) xmin : real(8), intent(in) xmax : real(8), intent(in)

—- 初期化 —-

[Source]

```  subroutine at_Initial(i,k,xmin,xmax)
integer, intent(in) :: i, k
real(8), intent(in) :: xmin,xmax

integer :: ii,kk

im=i
km=k
xl = xmax-xmin

if ( im <= 0 .or. km <= 0 ) then
call MessageNotify('E','at_initial',             'Number of grid points and waves should be positive')
elseif ( mod(im,2) /= 0 ) then
call MessageNotify('E','at_initial','Number of grid points should be even')
elseif ( km > im ) then
call MessageNotify('E','at_initial','KM shoud be less equal IM')
endif

allocate(t(3*im))
call fttcti(im,it,t)

allocate(g_X(0:im))
do ii=0,im
g_X(ii) = (xmax+xmin)/2 + xl/2 * cos(pi*ii/im)
enddo

allocate(g_X_Weight(0:im))
do ii=0,im
g_X_Weight(ii) = 1.0
do kk=2,km,2
g_X_Weight(ii) = g_X_Weight(ii)                           + 2/(1D0-kk**2) * cos(kk*ii*pi/im)
enddo
if ( (km == im) .and. (mod(im,2)==0) ) then  ! 最後の和は factor 1/2.
g_X_Weight(ii) = g_X_Weight(ii)                           - 1/(1D0-km**2)* cos(km*ii*pi/im)
endif
g_X_Weight(ii) = 2D0/im * g_X_Weight(ii) * xl/2
enddo
g_X_Weight(0)  = g_X_Weight(0) / 2
g_X_Weight(im) = g_X_Weight(im) / 2

end subroutine at_Initial
```