theta=linspace(0,2*pi,50); %theta variable
theta(length(theta))=[]; %remove the last one -- redundant
r=linspace(0,1,107); %r variable, odd number of values for integration later
a=cos(theta)'*r; % x-coordinate of points in domain
b=sin(theta)'*r; % y-coordinate of points in domain
%rho=0.3; %vertex of cone is at (rho,0)
denom1=-rho*(a-rho).^2+sqrt(rho^2*b.^4 ...
    +(b.^2+(a-rho).^2).*((a-rho).^2-rho^2*b.^2)); %with + square root
denom2=-rho*(a-rho).^2-sqrt(rho^2*b.^4 ...
    +(b.^2+(a-rho).^2).*((a-rho).^2-rho^2*b.^2)); % with - square root
f=(a-rho).*(b.^2+(a-rho).^2); %numerator of expression for cone
mask=(a>rho); % denominator switches based on sign of (a-rho)
f=1-f./(mask.*denom1+(~mask).*denom2); %choose right denominator
surf(a,b,f) %look at the cone
% columns of f are lines of constant r
% rows of f are lines of constant theta
fcoef=real(fft(f)); % rows are the Fourier coefficients as functions of r

fmax=5; %number of Fourier coefficients to expand in Bessel series
fcoeffapprox=zeros(fmax,length(r)); %to start building up approximate coefficients

%get Bessel series for first fmax Fourier coefficients
zinit=[2 5 8 11 14; 4 7 10 13 16; 5 8 11 15 18; 6 10 13 16 19; 7.5 10.5 14.3 17.6 20.8];
for j=0:fmax-1 %J0 series for c0, J1 series for c1, etc.
    
    for i=1:fmax
    z=zinit(j+1,i); %start finding zeros of Bessel functions by Newton's method
    incr=1;
    while(abs(incr)>1e-8)
        incr=2*besselj(j,z)/(besselj(j+1,z)-besselj(j-1,z));
        z=z+incr;
    end
    zr(j+1,i)=z;  %save accurate zero in array zr
    cf=SimpsonsRule(fcoef(j+1,:).*besselj(j,z*r).*r,0,1)/SimpsonsRule(besselj(j,z*r).^2.*r,0,1);
                  %find coefficient for Bessels series
    if(j==0)
        c(j+1,i)=cf;
    else
        c(j+1,i)=2*cf; %save coefficients of cos(j theta) in array c
    end
    fcoeffapprox(j+1,:)=fcoeffapprox(j+1,:)+cf*besselj(j,z*r); %add to Bessel series
    end
end

fcoeffapprox=[fcoeffapprox; zeros(length(theta)-2*fmax+1,length(r)); fcoeffapprox(fmax:-1:2,:)];
fapprox=real(ifft(fcoeffapprox)); %pad approximate coefficients, and ifft
%surf(a,b,fapprox)  %look at approximate cone

fs = 32768;  % Sampling rate (Hz)
timescale = 2*pi /fs; % Radians per sample
duration=2.0; %in seconds
t = timescale * (0:(duration*fs)); % scaled time for sinusoids

freqs=50*zr; %frequencies
initamp=c; % naive assignment of amplitudes?!
initamp=initamp/sum(sum(abs(initamp))); %limit max amplitude to 1
Q=300; % Quality of oscillators
damp=freqs/Q; % oscillators all with same Q
sz=size(damp);
szmat=sz(1)*sz(2); %total number of terms contributing to sound wave
damprow=zeros(1,szmat); %prepare row for parameters of sound wave terms
damprow(:)=damp; %put damping coefficients into the row
freqrow=zeros(1,szmat); %similarly for frequencies
freqrow(:)=freqs; %put frequencies into a row
initamprow=zeros(1,szmat); %similarly for amplitudes
initamprow(:)=initamp; %put amplitudes into a row
waves=exp(-damprow'*t).*sin(freqrow'*t); %each frequency, with its damping
tablawave=initamprow*waves; %add them up with appropriate amplitudes
sound(tablawave,fs); % play