function rayleigh_phase,ssa,nnormalize=nnormalize

    DELTA = (1.-.035)/(1.+.5*.035)
    PHASE = DELTA*.75*(1.+COS(SSA*!const.dtor)^2) + 1.-DELTA
;    phase=phase*0+1 ;This sets the Rayleigh phase to always be 1
    if not keyword_set(nnormalize) then phase=phase/(4.0d*!dpi)
    return,phase
end