Created
December 23, 2018 02:32
-
-
Save komasaru/368cefebfeb77fd2fcdf6d774ec4321a to your computer and use it in GitHub Desktop.
Fortran 95 source code to expand Fourier's series.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!**************************************************** | |
! フーリエ級数展開 | |
! f(t) = -1 (-PI < t <= 0 ) | |
! 1 ( 0 < t <= PI) | |
! | |
! date name version | |
! 2018.12.14 mk-mode.com 1.00 新規作成 | |
! | |
! Copyright(C) 2018 mk-mode.com All Rights Reserved. | |
!**************************************************** | |
! | |
module const | |
! SP: 単精度(4), DP: 倍精度(8) | |
integer, parameter :: SP = kind(1.0) | |
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP)) | |
real(DP), parameter :: PI = 4.0_DP * atan(1.0_DP) | |
end module const | |
module fourier | |
use const | |
implicit none | |
private | |
public :: expand | |
contains | |
! フーリエ級数展開 | |
! * 結果出力も行う | |
! | |
! :param(in) integer(4) n: 計算項数 | |
subroutine expand(n) | |
implicit none | |
integer(SP), intent(in) :: n | |
integer(SP) :: t, i | |
real(DP) :: y | |
y = 0.0_DP | |
do t = int(-1000.0_DP * PI), int(1000.0_DP * PI) | |
do i = 1, n | |
y = y + calc_term(i, t / 1.0e3_DP) | |
end do | |
print '(F6.3, X, F9.6)', t / 1.0e3_DP, 4.0_DP / PI * y | |
y = 0.0_DP | |
end do | |
end subroutine expand | |
! 各項計算 | |
function calc_term(n, t) result(c) | |
implicit none | |
integer(SP), intent(in) :: n | |
real(DP), intent(in) :: t | |
real(DP) :: c | |
c = sin((2 * n - 1) * t) / (2 * n - 1) | |
end function calc_term | |
end module fourier | |
program fourier_series_expansion | |
use const | |
use fourier | |
implicit none | |
integer(SP) :: n, t | |
! 計算項数の取得 | |
read (*, *) n | |
! フーリエ級数展開 | |
call expand(n) | |
end program fourier_series_expansion |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment