%
% COMSOL / Matlab Script
% Calculation of the acoustic radiation impedance of a flanged pipe 
% with the FEM (Helmholtz equation)
%
% Antoine Lefebvre, McGill University
%
% Copyright (c) 2010 Antoine Lefebvre
%

flclear fem

rho = 1.25;
c = 343;
a = 0.01;
L = 0.03;
N = 3; % ratio of sphere radius to cylinder radius

radiation_domain = ...
  block3(2*N*a,2*N*a,N*a,'base','corner','pos',[-N*a,-N*a,0]) * sphere3(N*a);
model = (cylinder3(a,L,'pos',[0,0,-L]) + radiation_somain) * ...
  block3(2*N*a,2*N*a,2*N*a,'pos',[0,0,-N*a]);

ka = [0.01 colon(0.1,0.1,1.0)];
k = ka/a;

clear s
s.objs={model};
s.name={'CO1'};
s.tags={'model'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Constants
fem.const = {'a','0.01 [m]', 'L','0.03 [m]', 'rho','1.25 [kg/m^3]', ...
             'c','343 [m/s]', 'v0','1 [m/s]', 'Zc','rho*c', 'S','pi*a^2'};

% Initialize mesh
edge = 0.001;
fem.mesh=meshinit(fem, 'hmax',0.005, 'hmaxfact',0.8, 'hcurve',0.5,  ...
                       'hgrad',1.3, 'hcutoff',0.02, 'hnarrow',0.6, ...
                       'hmaxedg',[12,edge]);

% Application mode 1
clear appl
appl.mode.class = 'AcoPressure';
appl.module = 'ACO';
appl.sshape = 3;
appl.assignsuffix = '_acpr';
clear prop
prop.elemdefault='Lag3';
appl.prop=prop;
clear bnd
bnd.p0 =  {0,0,0,1,0};
bnd.type = {'SH','cont','RAD','RAD','NA'};
bnd.wavetype = {'PL','PL','SPH','PL','PL'};
bnd.ind = [5,5,4,5,5,2,3,1,1];
appl.bnd = bnd;
clear equ
equ.rho = 'rho';
equ.cs = 'c';
equ.ind = [1,1];
appl.equ = equ;
appl.var = {'freq','freq'};
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

fem.expr = {'a0','i*omega_acpr*v0', 'freq','ka*c/(2*pi*a)'};

clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

fem=multiphysics(fem);
fem.xmesh=meshextend(fem);
fem.sol=femstatic(fem, 'solcomp',{'p'},'outcomp',{'p'}, 'pname','ka', ...
                       'plist', ka, 'oldcomp',{});
pin   =  postint(fem,'p/S',            'Dl',[3],'Edim',[2],...
                     'Solnum',[1:length(ka)]);
ZcUin = -postint(fem,'rho*c*nv_acpr/S','Dl',[3],'Edim',[2],...
                     'Solnum',[1:length(ka)]);

Zin = pin./ZcUin;
Zr = (Zin.*cos(k*L)-i*sin(k*L))./(-i*Zin.*sin(k*L)+cos(k*L));
flsave('flanged_pipe_quarter.mph', fem);
save('flanged_pipe_quarter.mat', 'Zr','a','ka');