Неклассифицированный оператор в (1) в математическом выражении

Question

Мой первый урок Fortran - построить функцию плотности вероятности радиальных штурмианских функций. В случае, если вас это интересует, радиальные штурмианские функции используются для построения собственных функций импульса для атома водорода.

Для того, чтобы произвести эти радиальные функции, необходимо сначала произвести некоторые многочлены, называемые полиномы Гегенбауэра, обозначается

С _б (х),

, где а и б должно быть сложены друг на друга. Нужно эти полиномы, потому что Sturmians (назовем их R_n, л) определяются как так,

R_n, л (р) = N ^{р л} ⁄ _{(р + к ) ^{л + 2}} С _{п - л - 1}^{л + 1} (^{р 2 - к} ⁄ _{р + к}),

, где N является константа нормировки, р-импульс, п является принцип квантового числа, л угловой момент и к является константой. Константа нормировки есть, так что, когда я прихожу к квадрату этой функции, она создаст распределение вероятности для импульса электрона в атоме водорода.

Гегенбауэра полиномы генерируются с помощью следующего рекуррентного соотношения:

C _п^л (х) = ⁄ _п [2 (л + п - 1) х С _{п - 1}^л (х) - (2l + п - 2) С _{п - 2}^л (х)],

с C ^л (х) = 1 и С 1 ^л (х) = 2LX, как вы могли заметить, л фиксируется, но п не является. В начале моей программы я укажу как l, так и n и выработаю полином Гегенбауэра, который мне нужен для радиальной функции, которую я хочу построить.

Проблемы, которые я имею с моим кодом на данный момент все в моей подпрограммой для разработки значение Гегенбауера полинома С _п-1^{л + 1} (^{р - к} ⁄ _{р + к}) для дополнительных значений р от 0 до 3.Я продолжаю получать ошибку

Unclassified statement at (1) 

но я не вижу, в чем проблема.

program Radial_Plot 

    implicit none 

    real, parameter :: pi = 4*atan(1.0) 
    integer, parameter :: top = 1000, l = 50, n = 100 
    real, dimension(1:top) :: x, y 
    real increment 
    real :: a=0.0, b = 2.5, k = 0.3 
    integer :: i 
    real, dimension(1:top) :: C 

    increment = (b-a)/(real(top)-1) 

    x(1) = 0.0 
    do i = 2, top 
     x(i) = x(i-1) + increment 
    end do 

    Call Gegenbauer(top, n, l, k, C) 

    y = x*C 
    ! y is the function that I shall be plotting between values a and b. 


end program Radial_Plot 


Subroutine Gegenbauer(top1, n1, l1, k1, CSub) 

    ! This subroutine is my attempt to calculate the Gegenbauer polynomials evaluated at a certain number of values between c and d. 

    implicit none 

    integer :: top1, i, j, n1, l1 
    real :: k1, increment1, c, d 
    real, dimension(1:top1) :: x1 
    real, dimension(1:n1 - l1, 1:top1) :: C1 
    real, dimension(1:n1 - l1) :: CSub 


    c = 0.0 
    d = 3.0 
    k1 = 0.3 
    n1 = 50 
    l1 = 25 
    top1 = 1000 


    increment1 = (d - c)/(real(top1) - 1) 

    x1(1) = 0.0 
    do i = 2, top1 
     x1(i) = x1(i-1) + increment1 
    end do 


    do j = 1, top1 

    C1(1,j) = 1 
    C1(2,j) = 2(l1 + 1)(x1(i)^2 - k1^2)/(x1(i)^2 + k1^2) 
    ! All the errors occurring here are all due to, and I quote, 'Unclassifiable statement at (1)', I can't see what the heck I have done wrong. 
     do i = 3, n1 - l1 
       C1(i,j) = 2(((l1 + 1)/n1) + 1)(x1(i)^2 - k1^2)/(x1(i)^2 + k1^2)C1(i,j-1) - ((2(l1+1)/n1) + 1)C1(i,j-2) 
     end do 

    CSub(j) = Cn(n1 - l1,j)^2 

    end do 
    return 
end Subroutine Gegenbauer 

Answer 1

Как francesalus правильно указал, что проблема в том, что вы используете ^ вместо ** для потенцирования. Кроме того, вы не ставите * между терминами, которые вы умножаете.

C1(1,j) = 1 
C1(2,j) = 2*(l1 + 1)*(x1(i)**2 - k1**2)/(x1(i)**2 + k1**2) 

do i = 3, n1 - l1 
      C1(i,j) = 2 * (((l1 + 1)/n1) + 1) * (x1(i)**2 - k1**2)/& 
         (x1(i)**2 + k1**2)*C1(i,j-1) - ((2(l1+1)/n1) + 1) * & 
         C1(i,j-2) 
end do 

CSub(j) = Cn(n1 - l1,j)**2 

С тех пор, как вы начинаете, у меня есть совет. Научитесь помещать все подпрограммы и функции в модули (если они не являются внутренними). Нет причины для оператора return в и подпрограмме, так же как заявление stop не требуется в программе и в программе.