bremsstrahlung Subroutine

subroutine bremsstrahlung
    !+ Calculates bremsstrahlung
    type(LocalProfiles) :: plasma
    integer :: i, ichan, nc, ic
    real(Float64) :: dlength, dlambda, gaunt, max_length
    real(Float64) :: spot_size, theta, sqrt_rho
    real(Float64), dimension(2) :: randomu
    real(Float64), dimension(3) :: vi, xyz, r0
    real(Float64), dimension(3,3) :: basis
    real(Float64), dimension(:), allocatable :: lambda_arr,brems

    allocate(lambda_arr(inputs%nlambda))
    allocate(brems(inputs%nlambda))

    do i=1,inputs%nlambda
        lambda_arr(i)= 10*((i-0.5)*inputs%dlambda+inputs%lambdamin) ! [A]
    enddo
    dlambda = 10*inputs%dlambda ![A]

    dlength = 0.3 !cm
    !! $OMP PARALLEL DO schedule(guided) private(ichan,xyz,vi,basis,spot_size, &
    !! $OMP& ic, nc,randomu,sqrt_rho,theta,r0,plasma,gaunt,brems)
    loop_over_channels: do ichan=1,spec_chords%nchan
        xyz = spec_chords%los(ichan)%lens
        vi = spec_chords%los(ichan)%axis
        vi = vi/norm2(vi)
        spot_size = spec_chords%los(ichan)%spot_size
        call line_basis(xyz,vi,basis)

        if(spot_size.le.0.d0) then
            nc = 1
        else
            nc = 100
        endif

        do ic=1,nc
            call randu(randomu)
            sqrt_rho = sqrt(randomu(1))
            theta = 2*pi*randomu(2)
            r0(1) = 0.d0
            r0(2) = spot_size*sqrt_rho*cos(theta)
            r0(3) = spot_size*sqrt_rho*sin(theta)
            r0 = matmul(basis,r0) + xyz

            ! Find edge of plasma
            call get_plasma(plasma,pos=r0)
            max_length=0.0
            do while (.not.plasma%in_plasma)
                r0 = r0 + vi*dlength ! move dlength
                call get_plasma(plasma,pos=r0)
                max_length = max_length + dlength
                if(max_length.gt.300) cycle loop_over_channels
            enddo

            ! Calculate bremsstrahlung along los
            do while (plasma%in_plasma)
                if(plasma%te.gt.0.0) then
                    gaunt = 5.542-(3.108-log(plasma%te))*(0.6905-0.1323/plasma%zeff)
                    brems = 7.57d-9*gaunt*plasma%dene**2*plasma%zeff/(lambda_arr &
                            *sqrt(plasma%te*1000.0))*exp(-h_planck*c0/(lambda_arr*plasma%te*1000.0)) &
                            *dlambda*(4.d0*pi)*1.d-4

                    spec%brems(:,ichan)= spec%brems(:,ichan) + (brems*dlength*1.d-2)/nc
                endif
                ! Take a step
                r0 = r0 + vi*dlength
                call get_plasma(plasma,pos=r0)
            enddo
        enddo

        if (inputs%verbose.ge.2)then
            WRITE(*,'(f7.2,"% completed",a,$)') 100*ichan/real(spec_chords%nchan),char(13)
        endif
    enddo loop_over_channels
    !! $OMP END PARALLEL DO

    deallocate(lambda_arr,brems)

end subroutine bremsstrahlung