pro dulk_gysy,delta,bb,theta,nv,freq,fi,rc,omega,tau


freq0=freq*1.e9

; Universal constant

c = 3e10         ; [cm/s] speed of light
kb = 1.38*10.^(-16)  ; [cgs] Bolzmann constant

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Dulk's approxmation limit
if ((theta lt 20) or (theta gt 90)) then begin
  print,'% Warning : assumed theta out of range for the approxmation'
  print,'%    It should be between  20 and 90 degree'
endif
if ((delta lt 2) or (delta gt 7)) then begin
  print,'% Warning : assumed delta out of range for the approxmation'
  print,'%    It should be between  2 and 7 '
endif

; calculation of e- cyclotron frequency eq(2)

fgyro = 2.8e6 *bb
fratio = freq0/fgyro

whr=where((fratio gt 100) or (fratio lt 10),count)
if (count ge 1) then begin
; Dulk's approxmation limit
  print,'% Warning : assumed B strength too weak or strong'
; print,'%    It should be between ',freq0/100./2.8e6, $
;   ' and ',freq0/10./2.8e6,' G'
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

sint = sin(theta*!dtor)
cost = cos(theta*!dtor)

; calculation of effective temperature  eq(37)

teeff = 2.2e9 * 10.^(-0.31*delta) * (sint)^(-0.36-0.06*delta) $
         * fratio^(0.50+0.085*delta)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; calculation of polarization eq(38)

rc = 1.26 * 10.^(0.035*delta-0.071*cost) *fratio^(-0.782+0.545*cost)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; calculation of flux

cm2_per_omega = (7.e7/(!dtor/3600.))^2
nvomega = nv / cm2_per_omega

; calculation of absorption coeff. eq(36)  
; calculation of optical thickness [tau = kappa *ll]
cgs2sfu = 1.e19 ; SFU/[erg/s/cm^2/Hz]

tauomega = 1.4*10.^(-9) * 10.^(-0.22*delta) * (sint)^(-0.09+0.72*delta) $
         * fratio^(-1.3-0.98*delta)  /bb *nvomega
fi = 2 * teeff * tauomega * kb *freq0^2 /c^2 * cgs2sfu


if keyword_set(omega) then tau = tauomega/omega


return
end