gauss_plume_second_order_rotated_func Function

    function gauss_plume_second_order_rotated_func(r,z,ay,by,az,bz,sig_y_0,sig_z_0,h)

        implicit none
        real r,z,ay,by,az,bz,sig_y_0,sig_z_0,h
        real gauss_plume_second_order_rotated_func
        real sig_th,sig_z,B,c
        real order_1,order_2
        real :: correction=2.
        real pi
        parameter (pi=3.141592)

        !Corrected for the B**2 falut in the taylor expansion and for the fact that the integral was only half a circle. 20.08.2019
        r=max(0.001,r)
        order_1=1.
        order_2=1.

        sig_th=(sig_y_0+ay*(exp(by*log(r))))/r
        sig_z=sig_z_0+az*(exp(bz*log(r)))
        !write(*,*) sig_z,sig_th*r,sig_z_0,sig_y_0
        B=-(sig_th**2)*(bz*(sig_z-sig_z_0)/r/sig_th+by*(r*sig_th-sig_y_0)/sig_z)
        !write(*,*) B
        if (B.gt.-1.) then
            !c=1./(2.*pi*sqrt(2.*pi)*r*sig_z*sqrt(1.+B))*tanh(2/sqrt(pi)*pi/(2.*sqrt(2.))/sig_th*sqrt(1.+B))*(exp((-(z-h)**2)/2./sig_z**2)+exp((-(z+h)**2)/2./sig_z**2))
            c=1./(2.*pi*sqrt(2.*pi)*r*sig_z*sqrt(1.+B))*erf(pi/(2.*sqrt(2.))/sig_th*sqrt(1.+B)*correction)*(exp((-(z-h)**2)/2./sig_z**2)+exp((-(z+h)**2)/2./sig_z**2))
        else
            c=correction/(4.*pi*sig_th*r*sig_z)*(1-order_1*pi**2*(1.+B)*correction**2/(24*sig_th**2)+order_2*pi**4*((1.+B)**2*correction**4/(640.*sig_th**4)))*(exp((-(z-h)**2)/2./sig_z**2)+exp((-(z+h)**2)/2./sig_z**2))
        endif

        gauss_plume_second_order_rotated_func=c

    end function gauss_plume_second_order_rotated_func