Waffle/src/WaveSystem.hpp at main · david-gaspard/Waffle · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
/****
 * @date Created on 2025-07-31 at 14:00:37 CEST
 * @author David Gaspard (ORCID 0000-0002-4449-8782) <david.gaspard@espci.fr>
 * @copyright This program is distributed under the MIT License.
 * @file C++ header providing the WaveSystem object, which describes a stationary scalar wave field defined on a SquareMesh.
 ***/
#ifndef _WAVE_SYSTEM_H
#define _WAVE_SYSTEM_H
#include "SquareMesh.hpp"
#include "SparseComplexMatrix.hpp"
#include "RealMatrix.hpp"

/**
 * Class defining the Square Mesh object.
 */
class WaveSystem {

    private:

    const std::string sysname;  // String containing the name of the system (typically describing the system geometry) which is used to generate file output.

    const SquareMesh mesh;   // Square mesh object.

    double kh;         // Wavenumber times the lattice step, 2*pi*h/lambda. Also the phase accumulated across a lattice step (in radian).
    double density;    // Density of scattering pixels per unit pixel. Value between 0 and 1.
    double holscat;    // Lattice step divided by the scattering mean free path, h/lscat. Total length is not a well defined unit.
    double holabso;    // Lattice step divided by the absorption length, h/labso. Total length is not a well defined unit.
    dcomplex khc;      // Complex wavenumber times the lattice step, khc = kh + I*(h/labso)/2.

    int npoint;        // Total number of points in the "mesh".
    int ninput;        // Number of points in the input, also equal to the number of input modes.
    int noutput;       // Number of points in the output, also equal to the number of output modes.
    int ninputprop;    // Number of propagating modes in the input lead(s).
    int noutputprop;   // Number of propagating modes in the output lead(s).
    double dosinput;   // Density of states in the input lead(s).
    double dosoutput;  // Density of states in the input lead(s).
    double doslattice; // Free density of states on a square lattice. Exact value based on the elliptic integral K(k).

    bool computed;     // Flag indicating if the Green function have been computed. "true" if "green" has been computed, "false" otherwise.
                       // This flag is set to "false" each time the Hamiltonian is modified by setDisorder(), and set to "true" by computeGreenFunction().

    // Internal matrices:
    SparseComplexMatrix hamiltonian;  // Hamiltonian matrix, i.e., discretization of the operator (d_x^2 + d_y^2 + k^2 + i*eps - U(x,y))*h^2
                                      // involved in the calculation of the Green function. Size: (npoint, npoint).
    SparseComplexMatrix inputState;   // Matrix storing the input modes in columns, and used as the independent term of the linear system
                                      // for the Green function. Size: (npoint, ninputprop).
    SparseComplexMatrix outputState;  // Matrix storing the output modes in columns. Size: (npoint, noutputprop).

    ComplexMatrix inputKlh;   // Column matrix containing the longitudinal wavenumbers in each channel. Size: (ninputprop, 1).
    ComplexMatrix outputKlh;  // Column matrix containing the longitudinal wavenumbers in each channel. Size: (noutputprop, 1).
    ComplexMatrix green;  // Matrix containing in columns the Green functions associated to each input modes. Size: (npoint, ninputprop).

    public:

    // Constructors/Destructors:
    WaveSystem(const std::string& name, const SquareMesh& mesh, const double kh, const double density, const double holscat, const double holabso);

    // Getters:
    int getNPoint() const;
    MeshPoint getPoint(const int ipoint) const;
    std::vector<Opening> getOpening() const;
    int getNOpening() const;
    int getNInput() const;
    int getNOutput() const;
    int getNInputProp() const;
    int getNOutputProp() const;
    int getNProp(const std::vector<int>& index) const;
    double getWavenumber() const;
    double getDensity() const;
    double getScattering() const;
    double getAbsorption() const;
    std::string getName() const;
    std::string uniqueFile(const std::string& dataname, const std::string& extension) const;
    std::vector<std::string> summary() const;

    // Print/save methods:
    void printSummary() const;
    void infoHamiltonian() const;
    void plotMatrixHamiltonian() const;
    void plotMatrixInputState() const;
    void plotMatrixOutputState() const;
    void plotMesh() const;
    void plotIntensity(const RealMatrix& intensity, const std::string& description, const std::string& dataname) const;
    void plotFields(const RealMatrix& fields, const std::vector<std::string>& labels, const std::string& description, const std::string& plotname) const;
    void plotGreenFunction();
    void plotInputModes(const int nmode);
    void plotTransmissionStates(const int nstate);

    // Public computational methods:
    void setDisorder(const uint64_t seed);
    void addPotential(const int x, const int y, const dcomplex uh2);
    void transmissionMatrix(ComplexMatrix& tmat);
    void reflectionMatrix(ComplexMatrix& rmat);
    void SBlockMatrix(const std::vector<int>& index, ComplexMatrix& smat);
    void checkUnitarity(const bool showtval);
    void checkResidual();
    void addTSpectrum(RealMatrix& tval);
    void addITransmission(const RealMatrix& trange, RealMatrix& tprofile, RealMatrix& nsample, RealMatrix& tval);
    void addITmax(RealMatrix& itmax, double& tmax);
    void addIIsotropic(RealMatrix& iavg);
    void addIMode(RealMatrix& iavg, const int imode);
    void addIPlane_v1(RealMatrix& iplane, double& tplane);
    void addIPlane_v2(RealMatrix& iplane, double& tplane);
    void addJTransmission(const RealMatrix& trange, RealMatrix& tcurrent, RealMatrix& nsample);
    void addPsiTransmission(const RealMatrix& trange, RealMatrix& tpsi, RealMatrix& nsample);
    void addDevTransmission(const RealMatrix& trange, RealMatrix& tdev, RealMatrix& nsample);
    void transmissionBalance(RealMatrix& tavg);

    private:

    // Private setters (because the Hamiltonian should be recomputed afterwards in principle):
    double checkWavenumber(const double kh) const;
    double checkDensity(const double density) const;
    double checkScattering(const double holscat) const;
    double checkAbsorption(const double holabso) const;
    dcomplex complexWavenumber(const double kh, const double holabso) const;

    // Private computational methods:
    void computeHamiltonian();
    int computeNInputProp() const;
    int computeNOutputProp() const;
    void computeIOStates();
    void wavefrontState(const std::vector<int>& index, SparseComplexMatrix& state, ComplexMatrix& klh, double& dos) const;
    void computeGreenFunction();
    void computeTransmissionStates(ComplexMatrix& tstate, RealMatrix& tval);
    void currentAt(const ComplexMatrix& psi, const int istate, const int ipoint, double& jx, double& jy) const;
    void deviationAt(const ComplexMatrix& psi, const int istate, const int ipoint, double& devx, double& devy, double& intensity) const;

};

#endif