program main

    external one,two,three,four

    integer :: n = 100
    real :: s1=0,s2=0,t1=0,t2=0,t3=0,t4=0,a=1,b=10 !ti - значения каждого интеграла
    real :: h=0
    real :: r(100)=0,theta(100)=0
    real :: phi(100,100)=0
    h=(b-a)/n
    s1=2./3. !cos(x)**2*sin(x) от нуля до пи
    s2=-8./15. !(1-3*cos(x)**2)*cos(x)**2*sin(x) от нуля до пи
    theta(1)=0
    theta(2)=acos(-1.0)/6.0
    theta(3)=acos(-1.0)/4.0
    theta(4)=acos(-1.0)/3.0
    theta(5)=acos(-1.0)/2.0
    theta(6)=acos(-1.0)

    open(11,file='phi.txt')
    do i=1,n !сетка
        r(i)=i*h
    end do

    do j=1,6 !для каждого угла
        do i=1,n
            call QG10(0,r(i),one,t1)
            call QG10(0, 1/r(i), two, t2)
            !call GL5(r(i),two,t2)
            !call GL5(r(i),four,t4)
            call QG10(0,r(i),three,t3)
            call QG10(0, 1/r(i), four, t4)
            phi(i,j)=2*acos(-1.0)*(s1*(t1/r(i)+t2)+(1-3*cos(theta(j))**2)*s2*(t3/r(i)**3+t4*r(i)**2)/4)
        end do
    end do
    do i=1,n
        write(11,*) r(i),phi(i,1),phi(i,2),phi(i,3),phi(i,4),phi(i,5),phi(i,6)
    end do
    write(*,*) t4
    close(11)
end

real function coeff(x)
    coeff=14440.0*x/(x**6-36*x**5+450*x**4-2400*x**3+5400*x**2-4320*x+720)**2
end function

subroutine GL5(r,f,s)
    real::x1=0.263560,x2=1.41340,x3=3.59643,x4=7.08581,x5=12.6408
    s=exp(-r)*(f(x1+r)*coeff(x1)+f(x2+r)*coeff(x2)+f(x3+r)*coeff(x3)+f(x4+r)*coeff(x4)+f(x5+r)*coeff(x5))
end subroutine

real function one(r)
    one=exp(-r)*r**2
end function

real function two(r)
    two=exp(-1./r)/r
end function

real function three(r)
    three=exp(-r)*r**4
end function

real function four(r)
    four=exp(-1./r)/r**3
end function

!real function two(r)
!    two=r
!end function

!real function four(r)
!    four=1.0/r
!end function