cylprop.f90

Go to the documentation of this file.

module cylprop
    use constants
    use solidprop
    use specfun
    implicit none
    private
    public trextcoeffintstrut,extcoeffintstrut,scattcoeffintstrut,cylconst
contains
!********************************BEGINNING*************************************
subroutine cylconst(lambda,theta,radius,Qs,Qt,gcyl)
    real(dp), intent(in) :: lambda
    real(dp), intent(in) :: theta
    real(dp), intent(in) :: radius
    real(dp), intent(out) :: Qs
    real(dp), intent(out) :: Qt
    real(dp), intent(out) :: gcyl
    integer :: &
        i,k,kmax,nwaw=1001
    real(dp) :: &
        x,&   !diffraction parameter
        tol=1e-10_dp,&   !tolerance for calculation of Bessel functions
        phi,& !scattered angle
        s1,s2,qte,qth,qse,qsh!,mue,muh
    complex(dp) :: &
        m,&    !complex index of refraction
        u,v,s,Dk,Rk,Ak1,Akm,Bk1,Bkm,Deltak
    complex(dp) :: &
        T1,T2,T3,T4
    real(dp), dimension(:), allocatable :: &
        adsr !angular distribution of scattered radiation
    complex(dp), dimension(:), allocatable :: &
        ae,ah,be,bh,ju,dju,yu,dyu,jv,djv,yv,dyv,h1v,dh1v,h2v,dh2v
    call optconst(lambda,n2,k2)
!    n2=1.6_dp
!    k2=0.01_dp
    m=cmplx(n2,-k2,kind=dp)
    x=2*pi*radius/lambda
    u=x*sqrt(m**2-sin(theta)**2)
    v=x*cos(theta)
    s=1/u**2-1/v**2
    kmax=nint(x+4*x**(1/3.0_dp)+2)+2 !+2 for safety (also there is most likely
    ! a mistake in Mackowski code on this place)
    if (kmax>50) kmax=50
!    write(*,*) kmax
    allocate(ae(0:kmax),ah(0:kmax),be(0:kmax),bh(0:kmax),ju(0:kmax),dju(0:kmax),&
        yu(0:kmax),dyu(0:kmax),jv(0:kmax),djv(0:kmax),yv(0:kmax),dyv(0:kmax),&
        h1v(0:kmax),dh1v(0:kmax),h2v(0:kmax),dh2v(0:kmax))
    call bessel(kmax,u,tol,ju,dju,yu,dyu)
    call bessel(kmax,v,tol,jv,djv,yv,dyv)
    call hankel(kmax,v,tol,h1v,dh1v,h2v,dh2v)
    do k=0,kmax
        dk=k*s*sin(theta)
        rk=jv(k)/h2v(k)
        ak1=dh2v(k)/(v*h2v(k))-dju(k)/(u*ju(k))
        akm=dh2v(k)/(v*h2v(k))-m**2*dju(k)/(u*ju(k))
        bk1=djv(k)/(v*jv(k))-dju(k)/(u*ju(k))
        bkm=djv(k)/(v*jv(k))-m**2*dju(k)/(u*ju(k))
        deltak=akm*ak1-dk**2
        ae(k)=iu*dk*rk*(bk1-ak1)/deltak
        ah(k)=rk*(akm*bk1-dk**2)/deltak
        be(k)=rk*(ak1*bkm-dk**2)/deltak
        bh(k)=-ae(k)
    enddo
    qte=real(0.5_dp*be(0),kind=dp)
    qth=real(0.5_dp*ah(0),kind=dp)
    qse=real(0.5_dp*be(0)*conjg(be(0)),kind=dp)
    qsh=real(0.5_dp*ah(0)*conjg(ah(0)),kind=dp)
    do k=1,kmax
        qte=qte+real(be(k),kind=dp)
        qth=qth+real(ah(k),kind=dp)
        qse=qse+real(ae(k)*conjg(ae(k))+be(k)*conjg(be(k)),kind=dp)
        qsh=qsh+real(ah(k)*conjg(ah(k))+bh(k)*conjg(bh(k)),kind=dp)
    enddo
    qte=4/x*qte
    qth=4/x*qth
    qse=4/x*qse
    qsh=4/x*qsh
    qs=(qse+qsh)/2
    qt=(qte+qth)/2
!    write(*,*) Qte,Qse,Qte-Qse
!    write(*,*) Qth,Qsh,Qth-Qsh
    !see Kaemmerlen, 2010: 10.1016/j.jqsrt.2009.11.018
!    Qs=Qs*(0.0015_dp*(lambda/(2*radius))**3-0.0219_dp*(lambda/(2*radius))**2+&
        ! 0.0688_dp*(lambda/(2*radius))+0.7639_dp)
!    Qt=Qt*(0.0015_dp*(lambda/(2*radius))**3-0.0219_dp*(lambda/(2*radius))**2+&
        ! 0.0688_dp*(lambda/(2*radius))+0.7639_dp)

    ! asymmetry factor is not needed, we need gcyl according to Placido,
    ! calculated as Dombrovsky
!    mue=Qse*sin(theta)**2*x/(4*cos(theta)**2)
!    muh=Qsh*sin(theta)**2*x/(4*cos(theta)**2)
!    do k=0,kmax-1
!        mue=mue+real(ae(k)*conjg(ae(k+1))+be(k)*conjg(be(k+1)),kind=dp)
!        muh=muh+real(ah(k)*conjg(ah(k+1))+bh(k)*conjg(bh(k+1)),kind=dp)
!    enddo
!    mue=mue*4*cos(theta)**2/x
!    muh=muh*4*cos(theta)**2/x
!    mu=(mue+muh)/(Qse+Qsh) !mue is actually mue*Qse

    allocate(adsr(nwaw))
    do i=1,nwaw
        phi=(i-1)*pi/(nwaw)
        t1=be(0)
        t2=ah(0)
        t3=0
        t4=0
        do k=1,kmax
            t1=t1+2*be(k)*cos(k*phi)
            t2=t2+2*ah(k)*cos(k*phi)
            t3=t3-2*iu*ae(k)*sin(k*phi)
            t4=-t3
        enddo
        adsr(i)=(abs(t1)**2+abs(t2)**2)/(pi*x*qs)
        adsr(i)=(abs(t1)**2+abs(t2)**2+abs(t3)**2+abs(t4)**2)/2
    enddo
    s1=adsr(1)*cos(0e0_dp)+adsr(nwaw)*cos(pi)
    s2=adsr(1)+adsr(nwaw)
    do i=2,nwaw-1
        phi=(i-1)*pi/(nwaw-1)
        s1=s1+2*adsr(i)*cos(phi)
        s2=s2+2*adsr(i)
    enddo
    s1=s1*pi/(2*(nwaw-1))
    s2=s2*pi/(2*(nwaw-1))
    gcyl=s1/s2
!    write(*,*) theta,s2*lambda/pi**2/radius,Qs
!    write(*,*) s1/s2,mu
!    phi=pi/2
!    T1=be(0)
!    T2=ah(0)
!    T3=0
!    T4=0
!    do k=1,kmax
!        T1=T1+2*be(k)*cos(k*phi)
!        T2=T2+2*ah(k)*cos(k*phi)
!        T3=T3-2*iu*ae(k)*sin(k*phi)
!        T4=-T3
!    enddo
!    write(*,*) abs(T1)**2/Qs,abs(T2)**2/Qs,abs(T3)**2/Qs

!    Qs=abs(be(0))**2+abs(ah(0))**2 !according to Modest, doesn't work
!    do i=1,kmax
!        Qs=Qs+abs(be(i))**2+abs(ah(i))**2+abs(bh(i))**2+abs(ae(i))**2
!    enddo
!    Qs=Qs/x
!    write(*,*) Qs
!    write(*,*) abs(T1)**2/Qs
    deallocate(adsr)
end subroutine cylconst
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function trextcoeffintstrut ( theta )
    real(dp), intent(in) :: theta
    real(dp) :: Qs,Qt,gcyl
    call cylconst(lambda,theta,dstrut/2,qs,qt,gcyl)
    trextcoeffintstrut=(qt-qs*sin(theta)**2-gcyl*qs*cos(theta)**2)*cos(theta)
end function trextcoeffintstrut
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function extcoeffintstrut ( theta )
    real(dp), intent(in) :: theta
    real(dp) :: Qs,Qt,gcyl
    call cylconst(lambda,theta,dstrut/2,qs,qt,gcyl)
    extcoeffintstrut=qt*cos(theta)
end function extcoeffintstrut
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function scattcoeffintstrut ( theta )
    real(dp), intent(in) :: theta
    real(dp) :: Qs,Qt,gcyl
    call cylconst(lambda,theta,dstrut/2,qs,qt,gcyl)
    scattcoeffintstrut=qs*cos(theta)
end function scattcoeffintstrut
!***********************************END****************************************
end module cylprop