/*
 * eigenmodes_in_muffler.java
 */

import com.comsol.model.*;
import com.comsol.model.util.*;

/** Model exported on May 15 2026, 10:09 by COMSOL 6.4.0.421. */
public class eigenmodes_in_muffler {

  public static Model run() {
    Model model = ModelUtil.create("Model");

//    From the File menu, choose New.
//    In the New window, click Model Wizard.
//    In the Model Wizard window, click 2D.
//    In the Select Physics tree, select Acoustics > Pressure Acoustics > Pressure Acoustics, Frequency Domain (acpr).
//    Click Add.
//    Click Study.
//    In the Select Study tree, select Preset Studies for Selected Physics Interfaces > Mode Analysis.
//    Click Done.

    model.component().create("comp1", true);

    model.component("comp1").geom().create("geom1", 2);

    model.component("comp1").mesh().create("mesh1");

    model.component("comp1").physics().create("acpr", "PressureAcoustics", "geom1");

    model.study().create("std1");
    model.study("std1").create("mode", "ModeAnalysis");
    model.study("std1").feature("mode").set("ftplistmethod", "manual");
    model.study("std1").feature("mode").set("chkeigregion", true);
    model.study("std1").feature("mode").set("storefact", false);
    model.study("std1").feature("mode").set("linpsolnum", "auto");
    model.study("std1").feature("mode").set("solnum", "auto");
    model.study("std1").feature("mode").set("notsolnum", "auto");
    model.study("std1").feature("mode").set("outputmap", new String[]{});
    model.study("std1").feature("mode").set("ngenAUX", "1");
    model.study("std1").feature("mode").set("goalngenAUX", "1");
    model.study("std1").feature("mode").set("ngenAUX", "1");
    model.study("std1").feature("mode").set("goalngenAUX", "1");

//    In the Model Builder window, under Component 1 (comp1), click Geometry 1.
//    In the Settings window for Geometry, locate the Units section.
//    From the Length unit list, select mm.

    model.component("comp1").geom("geom1").lengthUnit("mm");

//    In the Geometry toolbar, click Square.

    model.component("comp1").geom("geom1").create("sq1", "Square");

//    In the Settings window for Square, locate the Size section.
//    In the Side length text field, type 150.

    model.component("comp1").geom("geom1").feature("sq1").set("size", 150);

//    Locate the Position section.
//    From the Base list, select Center.

    model.component("comp1").geom("geom1").feature("sq1").set("base", "center");

//    In the Geometry toolbar, click Circle.

    model.component("comp1").geom("geom1").run("sq1");
    model.component("comp1").geom("geom1").create("c1", "Circle");

//    In the Settings window for Circle, locate the Size and Shape section.
//    In the Radius text field, type 75.

    model.component("comp1").geom("geom1").feature("c1").set("r", 75);

//    Locate the Position section.
//    In the x text field, type -75.

    model.component("comp1").geom("geom1").feature("c1").set("pos", new double[]{-75, 0});

//    In the Geometry toolbar, click Circle.

    model.component("comp1").geom("geom1").run("c1");
    model.component("comp1").geom("geom1").create("c2", "Circle");

//    In the Settings window for Circle, locate the Size and Shape section.
//    In the Radius text field, type 75.

    model.component("comp1").geom("geom1").feature("c2").set("r", 75);

//    Locate the Position section.
//    In the x text field, type 75.

    model.component("comp1").geom("geom1").feature("c2").set("pos", new double[]{75, 0});

//    In the Geometry toolbar, click Booleans and Partitions and choose Union.

    model.component("comp1").geom("geom1").run("c2");
    model.component("comp1").geom("geom1").create("uni1", "Union");

//    Click in the Graphics window and then press Ctrl+A to select all objects.

    model.component("comp1").geom("geom1").feature("uni1").selection("input").set("c1", "c2", "sq1");

//    In the Settings window for Union, locate the Union section.
//    Clear the Keep interior boundaries checkbox.

    model.component("comp1").geom("geom1").feature("uni1").set("intbnd", false);

//    Click Build All Objects.

    model.component("comp1").geom("geom1").runPre("fin");

//    Click the Zoom Extents button in the Graphics toolbar.

    model.component("comp1").geom("geom1").run();

//    In the Materials toolbar, click Add Material to open the Add Material window.
//    In the tree, select Built-in > Air.
//    Click Add to Component in the window toolbar.

    model.component("comp1").material().create("mat1", "Common");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("eta", "Piecewise");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("Cp", "Piecewise");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("rho", "Analytic");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("k", "Piecewise");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("cs", "Analytic");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("an1", "Analytic");
    model.component("comp1").material("mat1").propertyGroup("def").func().create("an2", "Analytic");
    model.component("comp1").material("mat1").propertyGroup()
         .create("RefractiveIndex", "RefractiveIndex", "Refractive index");
    model.component("comp1").material("mat1").propertyGroup()
         .create("NonlinearModel", "NonlinearModel", "Nonlinear model");
    model.component("comp1").material("mat1").propertyGroup().create("idealGas", "idealGas", "Ideal gas");
    model.component("comp1").material("mat1").propertyGroup("idealGas").func().create("Cp", "Piecewise");
    model.component("comp1").material("mat1").label("Air");
    model.component("comp1").material("mat1").set("family", "air");
    model.component("comp1").material("mat1").propertyGroup("def").func("eta").set("arg", "T");
    model.component("comp1").material("mat1").propertyGroup("def").func("eta")
         .set("pieces", new String[][]{{"200.0", "1600.0", "-8.38278E-7+8.35717342E-8*T^1-7.69429583E-11*T^2+4.6437266E-14*T^3-1.06585607E-17*T^4"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("eta").set("argunit", "K");
    model.component("comp1").material("mat1").propertyGroup("def").func("eta").set("fununit", "Pa*s");
    model.component("comp1").material("mat1").propertyGroup("def").func("Cp").set("arg", "T");
    model.component("comp1").material("mat1").propertyGroup("def").func("Cp")
         .set("pieces", new String[][]{{"200.0", "1600.0", "1047.63657-0.372589265*T^1+9.45304214E-4*T^2-6.02409443E-7*T^3+1.2858961E-10*T^4"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("Cp").set("argunit", "K");
    model.component("comp1").material("mat1").propertyGroup("def").func("Cp").set("fununit", "J/(kg*K)");
    model.component("comp1").material("mat1").propertyGroup("def").func("rho")
         .set("expr", "pA*0.02897/R_const[K*mol/J]/T");
    model.component("comp1").material("mat1").propertyGroup("def").func("rho").set("args", new String[]{"pA", "T"});
    model.component("comp1").material("mat1").propertyGroup("def").func("rho").set("fununit", "kg/m^3");
    model.component("comp1").material("mat1").propertyGroup("def").func("rho")
         .set("argunit", new String[]{"Pa", "K"});
    model.component("comp1").material("mat1").propertyGroup("def").func("rho")
         .set("plotaxis", new String[]{"off", "on"});
    model.component("comp1").material("mat1").propertyGroup("def").func("rho")
         .set("plotfixedvalue", new String[]{"101325", "273.15"});
    model.component("comp1").material("mat1").propertyGroup("def").func("rho")
         .set("plotargs", new String[][]{{"pA", "101325", "101325"}, {"T", "273.15", "293.15"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("k").set("arg", "T");
    model.component("comp1").material("mat1").propertyGroup("def").func("k")
         .set("pieces", new String[][]{{"200.0", "1600.0", "-0.00227583562+1.15480022E-4*T^1-7.90252856E-8*T^2+4.11702505E-11*T^3-7.43864331E-15*T^4"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("k").set("argunit", "K");
    model.component("comp1").material("mat1").propertyGroup("def").func("k").set("fununit", "W/(m*K)");
    model.component("comp1").material("mat1").propertyGroup("def").func("cs")
         .set("expr", "sqrt(1.4*R_const[K*mol/J]/0.02897*T)");
    model.component("comp1").material("mat1").propertyGroup("def").func("cs").set("args", new String[]{"T"});
    model.component("comp1").material("mat1").propertyGroup("def").func("cs").set("fununit", "m/s");
    model.component("comp1").material("mat1").propertyGroup("def").func("cs").set("argunit", new String[]{"K"});
    model.component("comp1").material("mat1").propertyGroup("def").func("cs")
         .set("plotfixedvalue", new String[]{"273.15"});
    model.component("comp1").material("mat1").propertyGroup("def").func("cs")
         .set("plotargs", new String[][]{{"T", "273.15", "373.15"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("an1").set("funcname", "alpha_p");
    model.component("comp1").material("mat1").propertyGroup("def").func("an1")
         .set("expr", "-1/rho(pA,T)*d(rho(pA,T),T)");
    model.component("comp1").material("mat1").propertyGroup("def").func("an1").set("args", new String[]{"pA", "T"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an1").set("fununit", "1/K");
    model.component("comp1").material("mat1").propertyGroup("def").func("an1")
         .set("argunit", new String[]{"Pa", "K"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an1")
         .set("plotaxis", new String[]{"off", "on"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an1")
         .set("plotfixedvalue", new String[]{"101325", "273.15"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an1")
         .set("plotargs", new String[][]{{"pA", "101325", "101325"}, {"T", "273.15", "373.15"}});
    model.component("comp1").material("mat1").propertyGroup("def").func("an2").set("funcname", "muB");
    model.component("comp1").material("mat1").propertyGroup("def").func("an2").set("expr", "0.6*eta(T)");
    model.component("comp1").material("mat1").propertyGroup("def").func("an2").set("args", new String[]{"T"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an2").set("fununit", "Pa*s");
    model.component("comp1").material("mat1").propertyGroup("def").func("an2").set("argunit", new String[]{"K"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an2")
         .set("plotfixedvalue", new String[]{"200"});
    model.component("comp1").material("mat1").propertyGroup("def").func("an2")
         .set("plotargs", new String[][]{{"T", "200", "1600"}});
    model.component("comp1").material("mat1").propertyGroup("def").set("thermalexpansioncoefficient", "");
    model.component("comp1").material("mat1").propertyGroup("def").set("molarmass", "");
    model.component("comp1").material("mat1").propertyGroup("def").set("bulkviscosity", "");
    model.component("comp1").material("mat1").propertyGroup("def")
         .set("thermalexpansioncoefficient", new String[]{"alpha_p(pA,T)", "0", "0", "0", "alpha_p(pA,T)", "0", "0", "0", "alpha_p(pA,T)"});
    model.component("comp1").material("mat1").propertyGroup("def").set("molarmass", "0.02897[kg/mol]");
    model.component("comp1").material("mat1").propertyGroup("def").set("bulkviscosity", "muB(T)");
    model.component("comp1").material("mat1").propertyGroup("def")
         .set("relpermeability", new String[]{"1", "0", "0", "0", "1", "0", "0", "0", "1"});
    model.component("comp1").material("mat1").propertyGroup("def")
         .set("relpermittivity", new String[]{"1", "0", "0", "0", "1", "0", "0", "0", "1"});
    model.component("comp1").material("mat1").propertyGroup("def").set("dynamicviscosity", "eta(T)");
    model.component("comp1").material("mat1").propertyGroup("def").set("ratioofspecificheat", "1.4");
    model.component("comp1").material("mat1").propertyGroup("def")
         .set("electricconductivity", new String[]{"0[S/m]", "0", "0", "0", "0[S/m]", "0", "0", "0", "0[S/m]"});
    model.component("comp1").material("mat1").propertyGroup("def").set("heatcapacity", "Cp(T)");
    model.component("comp1").material("mat1").propertyGroup("def").set("density", "rho(pA,T)");
    model.component("comp1").material("mat1").propertyGroup("def")
         .set("thermalconductivity", new String[]{"k(T)", "0", "0", "0", "k(T)", "0", "0", "0", "k(T)"});
    model.component("comp1").material("mat1").propertyGroup("def").set("soundspeed", "cs(T)");
    model.component("comp1").material("mat1").propertyGroup("def").addInput("temperature");
    model.component("comp1").material("mat1").propertyGroup("def").addInput("pressure");
    model.component("comp1").material("mat1").propertyGroup("RefractiveIndex")
         .set("n", new String[]{"1", "0", "0", "0", "1", "0", "0", "0", "1"});
    model.component("comp1").material("mat1").propertyGroup("NonlinearModel").set("BA", "def.gamma-1");
    model.component("comp1").material("mat1").propertyGroup("idealGas").func("Cp").label("Piecewise 2");
    model.component("comp1").material("mat1").propertyGroup("idealGas").func("Cp").set("arg", "T");
    model.component("comp1").material("mat1").propertyGroup("idealGas").func("Cp")
         .set("pieces", new String[][]{{"200.0", "1600.0", "1047.63657-0.372589265*T^1+9.45304214E-4*T^2-6.02409443E-7*T^3+1.2858961E-10*T^4"}});
    model.component("comp1").material("mat1").propertyGroup("idealGas").func("Cp").set("argunit", "K");
    model.component("comp1").material("mat1").propertyGroup("idealGas").func("Cp").set("fununit", "J/(kg*K)");
    model.component("comp1").material("mat1").propertyGroup("idealGas").set("Rs", "R_const/Mn");
    model.component("comp1").material("mat1").propertyGroup("idealGas").set("heatcapacity", "Cp(T)");
    model.component("comp1").material("mat1").propertyGroup("idealGas").set("ratioofspecificheat", "1.4");
    model.component("comp1").material("mat1").propertyGroup("idealGas").set("molarmass", "0.02897[kg/mol]");
    model.component("comp1").material("mat1").propertyGroup("idealGas").addInput("temperature");
    model.component("comp1").material("mat1").propertyGroup("idealGas").addInput("pressure");
    model.component("comp1").material("mat1").materialType("nonSolid");

//    In the Materials toolbar, click Add Material to close the Add Material window.
//    By default, the boundaries of the geometry will be considered to be sound hard walls. No other physics settings are needed.
//    In this model, the mesh is set up manually. Add a free triangular mesh and then proceed to the solver settings.
//    In the Mesh toolbar, click Free Triangular.

    model.component("comp1").mesh("mesh1").create("ftri1", "FreeTri");

//    In the Model Builder window, under Study 1, click Step 1: Mode Analysis.
//    In the Settings window for Mode Analysis, locate the Study Settings section.
//    Select the Desired number of modes checkbox.

    model.study("std1").feature("mode").set("neigsactive", true);

//    In the associated text field, type 8.

    model.study("std1").feature("mode").set("neigs", 8);

//    Select the Search for modes around shift checkbox.

    model.study("std1").feature("mode").set("shiftactive", true);

//    In the associated text field, type 20.

    model.study("std1").feature("mode").set("shift", "20");

//    The free-space propagation mode has an out-of-plane wave number equal to <c>omega/c</c> <c>=</c> <c>27.5</c> <c>rad/m</c>. With these settings, the solver returns the 8 modes with propagation constants closest to 20 rad/m first in the list.
//    In the Mode analysis frequency text field, type 1500[Hz].

    model.study("std1").feature("mode").set("modeFreq", "1500[Hz]");

//    This setting makes the software look for propagating modes with cutoff frequencies up to 1500 Hz.
//    In the Study toolbar, click Compute.

    model.study("std1").createAutoSequences("all");

    model.sol("sol1").runAll();

    model.result().create("pg1", "PlotGroup2D");
    model.result("pg1").set("data", "dset1");
    model.result("pg1").setIndex("looplevel", 1, 0);
    model.result("pg1").create("surf1", "Surface");
    model.result("pg1").feature("surf1").set("expr", new String[]{"acpr.p_t"});
    model.result("pg1").feature("surf1").set("colortable", "Wave");
    model.result("pg1").feature("surf1").set("colorscalemode", "linearsymmetric");
    model.result("pg1").set("showlegendsunit", true);
    model.result("pg1").feature("surf1").set("colortable", "WaveLight");
    model.result("pg1").create("con1", "Contour");
    model.result("pg1").feature("con1").set("expr", new String[]{"acpr.p_t"});
    model.result("pg1").feature("con1").set("colortable", "Wave");
    model.result("pg1").feature("con1").set("colorscalemode", "linearsymmetric");
    model.result("pg1").feature("con1").set("colorlegend", false);
    model.result("pg1").label("Acoustic Pressure (acpr)");
    model.result().create("pg2", "PlotGroup2D");
    model.result("pg2").set("data", "dset1");
    model.result("pg2").setIndex("looplevel", 1, 0);
    model.result("pg2").create("surf1", "Surface");
    model.result("pg2").feature("surf1").set("expr", new String[]{"acpr.Lp_t"});
    model.result("pg2").feature("surf1").set("colortable", "Rainbow");
    model.result("pg2").feature("surf1").set("colorscalemode", "linear");
    model.result("pg2").set("showlegendsunit", true);
    model.result("pg2").label("Sound Pressure Level (acpr)");
    model.result("pg1").run();

//    The solver has found the free-space mode and all other propagating modes. There is a total of 5 different propagating modes. Because the waves can propagate both into and out of the modeling plane, each mode gets reported twice, with positive and negative out-of-plane wave numbers.
//    For the positive out-of-plane wave numbers, it holds that the higher the mode, the lower the wave number. However, the solver does not stop at zero. Because you asked for more than the 5 existing propagating modes, you get additional modes with imaginary out-of-plane wave numbers. This indicates that they are evanescent. The default plot shows the acoustic pressure distribution for a mode with a wave number of -17.95i rad/m.
//    In the Settings window for 2D Plot Group, locate the Data section.
//    From the Out-of-plane wave number (rad/m) list, select -15.461.

    model.result("pg1").set("looplevel", new int[]{4});

//    In the Acoustic Pressure (acpr) toolbar, click Plot.

    model.result("pg1").run();

//    Click the Zoom Extents button in the Graphics toolbar.
//    In the Results toolbar, click Evaluation Group.

    model.result().evaluationGroup().create("eg1", "EvaluationGroup");

//    In the Settings window for Evaluation Group, type Evaluation Group 1 - Mode Cutoff Frequency in the Label text field.

    model.result().evaluationGroup("eg1").label("Evaluation Group 1 - Mode Cutoff Frequency");

//    Locate the Data section.
//    From the Out-of-plane wave number selection list, select From list.

    model.result().evaluationGroup("eg1").setIndex("looplevelinput", "manual", 0);

//    In the Out-of-plane wave number (rad/m) list, select -15.461.

    model.result().evaluationGroup("eg1").setIndex("looplevel", new int[]{4}, 0);

//    select -15.461 in the Out-of-plane wave number (rad/m) list.
//    Right-click Evaluation Group 1 - Mode Cutoff Frequency and choose Point Evaluation.

    model.result().evaluationGroup("eg1").create("pev1", "EvalPoint");

//    Select Point 1.

    model.result().evaluationGroup("eg1").feature("pev1").selection().set(1);

//    In the Settings window for Point Evaluation, locate the Expressions section.
//    In the table, enter the following settings:

    model.result().evaluationGroup("eg1").feature("pev1")
         .setIndex("expr", "sqrt(acpr.omega^2-acpr.kz^2*acpr.c^2)/(2*pi)", 0);
    model.result().evaluationGroup("eg1").feature("pev1").setIndex("descr", "Mode Cutoff Frequency", 0);

//    In the Evaluation Group 1 - Mode Cutoff Frequency toolbar, click Evaluate.

    model.result().evaluationGroup("eg1").run();

//    Finally, create a plot to visualize the mode in 3D using the <l>2D Extrusion</l> dataset.
//    In the Results toolbar, click More Datasets and choose Extrusion 2D.

    model.result().dataset().create("extr1", "Extrude2D");

//    In the Settings window for Extrusion 2D, locate the Extrusion section.
//    In the z maximum text field, type 1000.

    model.result().dataset("extr1").set("zmax", "1000");

//    In the Resolution text field, type 100.

    model.result().dataset("extr1").set("res", 100);

//    Click to expand the Advanced section.
//    In the Out-of-plane wave number text field, type acpr.kz.

    model.result().dataset("extr1").set("wavenumber", "acpr.kz");

//    In the Results toolbar, click 3D Plot Group.

    model.result().create("pg3", "PlotGroup3D");
    model.result("pg3").run();

//    In the Settings window for 3D Plot Group, type Mode 3D in the Label text field.

    model.result("pg3").label("Mode 3D");

//    Locate the Data section.
//    From the Out-of-plane wave number (rad/m) list, select -15.461.

    model.result("pg3").set("looplevel", new int[]{4});

//    Right-click Mode 3D and choose Surface.

    model.result("pg3").create("surf1", "Surface");
    model.result("pg3").feature("surf1").set("evaluationsettings", "parent");
    model.result("pg3").run();

//    Rotate and zoom to get the desired 3D view.

    model.title("Eigenmodes in a Muffler");

    model
         .description("This example consists of a 2D analysis of propagation modes in the chamber of the muffler used in the Absorptive Muffler model.");

    return model;
  }

  public static void main(String[] args) {
    run();
  }

}