gauss_plume_second_order_rotated_integral_func Function

    function gauss_plume_second_order_rotated_integral_func(r,ay,by,az,bz,sig_y_0,sig_z_0,h,H1,H2)

        implicit none
        real r,ay,by,az,bz,sig_y_0,sig_z_0,h,H1,H2
        real gauss_plume_second_order_rotated_integral_func
        real sig_th,sig_z,B,c_y_int,c_z_int
        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
        !Still need to implement reflections
        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)))

        B=-(sig_th**2)*(bz*(sig_z-sig_z_0)/r/sig_th+by*(r*sig_th-sig_y_0)/sig_z)

        c_z_int=sqrt(pi/2.)*sig_z*(erf((H2-h)/sqrt(2.)/sig_z)-erf((H1-h)/sqrt(2.)/sig_z)+erf((H2+h)/sqrt(2.)/sig_z)-erf((H1+h)/sqrt(2.)/sig_z))/(H2-H1)

        if (B.gt.-1.) then
            !c_int=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))
            c_y_int=1./(2.*pi*sqrt(2.*pi)*r*sig_z*sqrt(1.+B))*erf(pi/(2.*sqrt(2.))/sig_th*sqrt(1.+B)*correction)
        else
            c_y_int=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)))
        endif

        gauss_plume_second_order_rotated_integral_func=c_y_int*c_z_int

    end function gauss_plume_second_order_rotated_integral_func