That’s the code of the UDO I have, working on macOS and on Windows (I developed my csd there):
/*
m_aerosound - Aeroacoustic semi-physical model
DESCRIPTION
Model the effect of wind blowing past a cylinder.
SYNTAX
aleft, aright m_aerosound kairspeed, kdiametre, klength, kelevation, \
kazimuth, kdistance
INITIALIZATION
PERFORMANCE
kairspeed - speed of the wind in m/s (may not be 0!)
kdiametre - diametre of the cylinder in m (may not be 0!)
klength - length of the cylinder in metres
kelevation - elevation angle between listener and sound source (in Radians)
kazimuth - Azimuth angle between listener and sound source (in Radians)
kdistance - Distance between listener and sound source (in Metres)
CREDITS
Written in Csound by Jeanette C., original proposition, research and PureData
code by Rod Selfridge as part of his thesis at Queen Mary University of
London, many thanks to him for his untiring support and generousity!
*/
opcode m_aerosound, aa, kkkkkk
kairspeed, kdiametre, klength, kelevation, kazimuth, kdistance xin
; Protect against division by 0, airspeed must be greater than 0
if kairspeed == 0 then
kairspeed = 1
endif
kbasefreq = 0.2 * kairspeed / kdiametre ; base frequency of the Aeolian tone
kM = kairspeed / 343 ; the Mach number
; Calculate the Reynolds number necessary as parameter to the q-value
kreynolds = 1.225 * kdiametre * kairspeed / 0.000148
; Calculate the inverse of the Q-value percentage, based on value of Reynolds
if kreynolds < 193260 then
kqinv = 0.00004624 * kreynolds + 0.9797
else
kqinv = 0.000000000127 * (kreynolds^2) - 0.00008522 * kreynolds + 16.5
endif
kq = 1 / (kqinv * 0.01) ; invert the percentage
kq port kq, 0.5 ; smooth the q-value, in case of passing the reynolds
; threshhold in the if-statement above
kbaseamp = 0.000000003 * (klength * (kairspeed^6) * (sin(kelevation)^2) \
* (cos(kazimuth)^2)) / ((kdistance^2) * ((1 - kM * cos(kelevation)^4)))
; amplitude of the base frequency, the first factor is a static
; correction, so we don't overflow and distort
kampfactor = log10(kbaseamp / 2 * 0.00001) ; linear/log conversion factor
; present in all amplitude calculation for the harmonics
; Amplitudes of the harmonics for drag and lift dipoles
k2ndamp = 10^(0.1 * kampfactor)
k3rdamp = 10^(0.6 * kampfactor)
k4thamp = 10^(0.0125 * kampfactor)
k5thamp = 10^(0.1 * kampfactor)
anoise fractalnoise 1, 1 ; Brown noise as basic sound source
; All noise bands for basic tone and harmonics
abase butterbp anoise*kbaseamp, kbasefreq, kq
a2nd butterbp anoise*k2ndamp, kbasefreq*2, kq
a3rd butterbp anoise*k3rdamp, kbasefreq*3, kq
a4th butterbp anoise*k4thamp, kbasefreq*4, kq
a5th butterbp anoise*k5thamp, kbasefreq*5, kq
aout = abase + a2nd + a3rd + a4th + a5th
aout butterhp aout, kbasefreq
; Convert angles to polar coordinates
kx = sin(kazimuth) * kdistance
ky = cos(kazimuth) * kdistance
kz = sin(kelevation) * kdistance
; Place the sound in 3d space with b-format output
aw, ax, ay, az spat3d aout, kx, ky, kz, 1, 0, 3, 1, 2
; Encode b-format to stereo
aout_l, aout_r bformdec1 1, aw, ax, ay, az
xout aout_l, aout_r
endop