************80
! FCOEF: expansion coefficients for Mathieu and modified Mathieu functions.
Licensing:
This routine is copyrighted by Shanjie Zhang and Jianming Jin. However,
they give permission to incorporate this routine into a user program
provided that the copyright is acknowledged.
Modified:
01 August 2012
Author:
Shanjie Zhang, Jianming Jin
Reference:
Shanjie Zhang, Jianming Jin,
Computation of Special Functions,
Wiley, 1996,
ISBN: 0-471-11963-6,
LC: QA351.C45.
Parameters:
Input, integer ( kind = 4 ) KD, the case code.
1, for cem(x,q) ( m = 0,2,4,...)
2, for cem(x,q) ( m = 1,3,5,...)
3, for sem(x,q) ( m = 1,3,5,...)
4, for sem(x,q) ( m = 2,4,6,...)
Input, integer ( kind = 4 ) M, the order of the Mathieu function.
Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.
Input, real ( kind = 8 ) A, the characteristic value of the Mathieu
functions for given m and q.
Output, real ( kind = 8 ) FC(*), the expansion coefficients of Mathieu
functions ( k = 1,2,...,KM ). FC(1),FC(2),FC(3),... correspond to
A0,A2,A4,... for KD = 1 case,
A1,A3,A5,... for KD = 2 case,
B1,B3,B5,... for KD = 3 case,
B2,B4,B6,... for KD = 4 case.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer(kind=4) | :: | kd | ||||
integer(kind=4) | :: | m | ||||
real(kind=8) | :: | q | ||||
real(kind=8) | :: | a | ||||
real(kind=8) | :: | fc(251) |
subroutine fcoef ( kd, m, q, a, fc ) !*****************************************************************************80 ! !! FCOEF: expansion coefficients for Mathieu and modified Mathieu functions. ! ! Licensing: ! ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, ! they give permission to incorporate this routine into a user program ! provided that the copyright is acknowledged. ! ! Modified: ! ! 01 August 2012 ! ! Author: ! ! Shanjie Zhang, Jianming Jin ! ! Reference: ! ! Shanjie Zhang, Jianming Jin, ! Computation of Special Functions, ! Wiley, 1996, ! ISBN: 0-471-11963-6, ! LC: QA351.C45. ! ! Parameters: ! ! Input, integer ( kind = 4 ) KD, the case code. ! 1, for cem(x,q) ( m = 0,2,4,...) ! 2, for cem(x,q) ( m = 1,3,5,...) ! 3, for sem(x,q) ( m = 1,3,5,...) ! 4, for sem(x,q) ( m = 2,4,6,...) ! ! Input, integer ( kind = 4 ) M, the order of the Mathieu function. ! ! Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions. ! ! Input, real ( kind = 8 ) A, the characteristic value of the Mathieu ! functions for given m and q. ! ! Output, real ( kind = 8 ) FC(*), the expansion coefficients of Mathieu ! functions ( k = 1,2,...,KM ). FC(1),FC(2),FC(3),... correspond to ! A0,A2,A4,... for KD = 1 case, ! A1,A3,A5,... for KD = 2 case, ! B1,B3,B5,... for KD = 3 case, ! B2,B4,B6,... for KD = 4 case. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) f real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) f3 real ( kind = 8 ) fc(251) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) kb integer ( kind = 4 ) kd integer ( kind = 4 ) km integer ( kind = 4 ) l integer ( kind = 4 ) m real ( kind = 8 ) q real ( kind = 8 ) qm real ( kind = 8 ) s real ( kind = 8 ) s0 real ( kind = 8 ) sp real ( kind = 8 ) ss real ( kind = 8 ) u real ( kind = 8 ) v if ( q <= 1.0D+00 ) then qm = 7.5D+00 + 56.1D+00 * sqrt ( q ) - 134.7D+00 * q & + 90.7D+00 * sqrt ( q ) * q else qm = 17.0D+00 + 3.1D+00 * sqrt ( q ) - 0.126D+00 * q & + 0.0037D+00 * sqrt ( q ) * q end if km = int ( qm + 0.5D+00 * m ) if ( q == 0.0D+00 ) then do k = 1, km fc(k) = 0.0D+00 end do if ( kd == 1 ) then fc((m+2)/2) = 1.0D+00 if (m == 0 ) then fc(1) = 1.0D+00 / sqrt ( 2.0D+00 ) end if else if ( kd == 4 ) then fc(m/2) = 1.0D+00 else fc((m+1)/2) = 1.0D+00 end if return end if kb = 0 s = 0.0D+00 f = 1.0D-100 u = 0.0D+00 fc(km) = 0.0D+00 if ( kd == 1 ) then l = 0 do k = km, 3, -1 v = u u = f f = ( a - 4.0D+00 * k * k ) * u / q - v if ( abs ( f ) < abs ( fc(k+1) ) ) then kb = k fc(1) = 1.0D-100 sp = 0.0D+00 f3 = fc(k+1) fc(2) = a / q * fc(1) fc(3) = ( a - 4.0D+00 ) * fc(2) / q - 2.0D+00 * fc(1) u = fc(2) f1 = fc(3) do i = 3, kb v = u u = f1 f1 = ( a - 4.0D+00 * ( i - 1.0D+00 ) ** 2 ) * u / q - v fc(i+1) = f1 if ( i == kb ) then f2 = f1 else sp = sp + f1 * f1 end if end do sp = sp + 2.0D+00 * fc(1) ** 2 + fc(2) ** 2 + fc(3) ** 2 ss = s + sp * ( f3 / f2 ) ** 2 s0 = sqrt ( 1.0D+00 / ss ) do j = 1, km if ( j <= kb + 1 ) then fc(j) = s0 * fc(j) * f3 / f2 else fc(j) = s0 * fc(j) end if end do l = 1 exit else fc(k) = f s = s + f * f end if end do if ( l == 0 ) then fc(2) = q * fc(3) / ( a - 4.0D+00 - 2.0D+00 * q * q / a ) fc(1) = q / a * fc(2) s = s + 2.0D+00 * fc(1) ** 2 + fc(2) ** 2 s0 = sqrt ( 1.0D+00 / s ) do k = 1, km fc(k) = s0 * fc(k) end do end if else if ( kd == 2 .or. kd == 3 ) then l = 0 do k = km, 3, -1 v = u u = f f = ( a - ( 2.0D+00 * k - 1 ) ** 2 ) * u / q - v if ( abs ( fc(k) ) <= abs ( f ) ) then fc(k-1) = f s = s + f * f else kb = k f3 = fc(k) l = 1 exit end if end do if ( l == 0 ) then fc(1) = q / ( a - 1.0D+00 - ( - 1 ) ** kd * q ) * fc(2) s = s + fc(1) * fc(1) s0 = sqrt ( 1.0D+00 / s ) do k = 1, km fc(k) = s0 * fc(k) end do else fc(1) = 1.0D-100 fc(2) = ( a - 1.0D+00 - ( - 1 ) ** kd * q ) / q * fc(1) sp = 0.0D+00 u = fc(1) f1 = fc(2) do i = 2, kb - 1 v = u u = f1 f1 = ( a - ( 2.0D+00 * i - 1.0D+00 ) ** 2 ) * u / q - v if ( i /= kb - 1 ) then fc(i+1) = f1 sp = sp + f1 * f1 else f2 = f1 end if end do sp = sp + fc(1) ** 2 + fc(2) ** 2 ss = s + sp * ( f3 / f2 ) ** 2 s0 = 1.0D+00 / sqrt ( ss ) do j = 1, km if ( j < kb ) then fc(j) = s0 * fc(j) * f3 / f2 else fc(j) = s0 * fc(j) end if end do end if else if ( kd == 4 ) then l = 0 do k = km, 3, -1 v = u u = f f = ( a - 4.0D+00 * k * k ) * u / q - v if ( abs ( fc(k) ) <= abs ( f ) ) then fc(k-1) = f s = s + f * f else kb = k f3 = fc(k) l = 1 exit end if end do if ( l == 0 ) then fc(1) = q / ( a - 4.0D+00 ) * fc(2) s = s + fc(1) * fc(1) s0 = sqrt ( 1.0D+00 / s ) do k = 1, km fc(k) = s0 * fc(k) end do else fc(1) = 1.0D-100 fc(2) = ( a - 4.0D+00 ) / q * fc(1) sp = 0.0D+00 u = fc(1) f1 = fc(2) do i = 2, kb - 1 v = u u = f1 f1 = ( a - 4.0D+00 * i * i ) * u / q - v if ( i /= kb - 1 ) then fc(i+1) = f1 sp = sp + f1 * f1 else f2 = f1 end if end do sp = sp + fc(1) ** 2 + fc(2) ** 2 ss = s + sp * ( f3 / f2 ) ** 2 s0 = 1.0D+00 / sqrt ( ss ) do j = 1, km if ( j < kb ) then fc(j) = s0 * fc(j) * f3 / f2 else fc(j) = s0 * fc(j) end if end do end if end if if ( fc(1) < 0.0D+00 ) then do j = 1, km fc(j) = -fc(j) end do end if return end subroutine fcoef