# Single Scattering Albedo The single scattering albedo is $$a=1$$ , and the total optical depth $$\tau^{\star}=0.01$$ . The setting satisfies single-scattering approximation, $$a \times \tau^{\star}=0.01\ll1$$. The incident solar flux is an unpolarized light, $$F=\[0.5, 0.5, 0, 0]$$. In the following case example, the angles are set as follows, $$\theta^{\prime}=0^{o}, \mu^{\prime}=1, \phi^{\prime}=0^{o} ; \theta=\left\[0^{o}, 180^{\circ}\right], \mu=\[-1,1], \phi=60^{o}$$, where $$\theta^{\prime}, \phi^{\prime}$$ are the polar angle and azimuth angle of the incident light $$\Omega^{\prime}$$ , and the $$\theta , \phi$$ are for the corresponding observe angles $$\Omega$$ . ```python a = 1; tau = 0.01; F = [.5;.5;0;0]; thetap = 0.00001; mup = cos(thetap*pi/180); phip = 0; phi = 60; ``` Compute the Legendre-Gauss nodes ( $$ctheta$$ ) and weights ( $$Weight$$ ) on an interval \[-1,1] with truncation order nthetas ( $$ntheta=400$$ ) using an outside function `lgwt`, and retrieve the corresponding $$thetas$$. Function `lgwt` takes integral using Legendre-Gauss Quadrature. ```python nthetas = 400; [ctheta, Weight] = lgwt(nthetas,-1,1); theta = acos(ctheta)*180/pi; ``` The angle $$\Theta$$ between the directions of incidence $$\Omega$$ and observation $$\Omega^{\prime}$$ is given by $$ \cos \Theta=\cos \theta^{\prime} \cos \theta+\sin \theta^{\prime} \sin \theta \cos \left(\phi^{\prime}-\phi\right) $$ ![Illustration of the relationship between Cartesian and spherical coordinates.](https://492391472-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LiBoWdc5sh0EFrazIOe%2F-LlhCUnL2WJX7FGUhVaN%2F-Llhq082K7WWtT3-pEvM%2FCartesian.png?alt=media\&token=dcc60b83-7ebb-4157-af2d-d897d8084681) There are two representations for the phase matrix, $$\mathbf{P}*{S}$$ and $$\mathbf{P}*{C}$$, where $$\mathbf{P}*{S}$$ is used for the Stocks vector representation $$\mathbf{I}*{S}=\[I, Q, U, V]^{T}$$, and $$\mathbf{P}*{C}$$ is for $$\mathbf{I}*{C}=\left\[I\_{ |}, I\_{\perp}, U, V\right]^{T}$$ representation (or Chandrasekhar's representation). The connection between these two representations is simply $$\mathbf{I}*{S} = \mathbf{DI}*{C}$$ , where the matrix $$\mathbf{D}$$ is given by: $$ \mathbf{D} \equiv\left(\begin{array}{cccc}{1} & {1} & {0} & {0} \ {1} & {-1} & {0} & {0} \ {0} & {0} & {1} & {0} \ {0} & {0} & {0} & {1}\end{array}\right) $$ The phase matrix $$\mathbf{P}*{C}$$ in $$\mathbf{I}*{C}=\left\[I\_{ |}, I\_{\perp}, U, V\right]^{T}$$ representation is related to the phase matrix $$\mathbf{P}*{S}$$ in $$\mathbf{I}*{S}=\[I, Q, U, V]^{T}$$ as $$\mathbf{P}*{C}=\mathbf{D}^{-1} \mathbf{P}*{S} \mathbf{D}$$, and $$\mathbf{D}^{-1}$$is given by: $$ \mathbf{D}^{-1}=\left(\begin{array}{cccc}{\frac{1}{2}} & {\frac{1}{2}} & {0} & {0} \ {\frac{1}{2}} & {-\frac{1}{2}} & {0} & {0} \ {0} & {0} & {1} & {0} \ {0} & {0} & {0} & {1}\end{array}\right) $$ ```python cTheta = mup*ctheta+sqrt((1-mup^2)*(1-ctheta.^2))*cos((phip-phi)/180*pi); Theta = acos(cTheta)*180/pi; D = [1 1 0 0;... 1 -1 0 0;... 0 0 1 0;... 0 0 0 1]; ``` In this test case, we only include one wavelength at 500 nm, number of moment 300+1, and load the six columns of greek constants from `Mie_tool output` file, without delta fit. ```python number of wavelength numwl = 1; % number of moment nmom = 300+1; fname1 = 'mie_output_small2.dat'; fid1 = fopen(fname1); rawdata1(1, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 3); % rawdata1(2, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); % rawdata1(3, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); % rawdata1(4, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); % rawdata1(5, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); % rawdata1(6, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); % rawdata1(7, :) = textscan(fid1, '%f %f %f %f %f %f', nmom, 'headerLines', 4); fclose(fid1); rawdata2 = cell2mat(rawdata1); µdata2 = zeros(nmom, 6, numwl);µ ``` The expression for the Fourier coefficients can be written as $$ \mathbf{A}^{m}\left(\tau, u^{\prime}, u\right)=\sum\_{\ell=m}^{2 N-1} \mathbf{P}*{\ell}^{m}(u) \boldsymbol{\Lambda}*{\ell}(\tau) \mathbf{P}\_{\ell}^{m}\left(u^{\prime}\right) $$ where $$ \boldsymbol{\Lambda}*{\ell}(\tau)=\left(\begin{array}{cccc}{\alpha*{1}^{\ell}} & {\beta\_{1}^{\ell}} & {0} & {0} \ {\beta\_{1}^{\ell}} & {\alpha\_{2}^{\ell}} & {0} & {0} \ {0} & {0} & {\alpha\_{3}^{\ell}} & {\beta\_{2}^{\ell}} \ {0} & {0} & {-\beta\_{2}^{\ell}} & {\alpha\_{4}^{\ell}}\end{array}\right) $$ The matrices $$P\_{\ell}^{m}( \pm \mu)$$ are defined through $$ \mathbf{P}*{\ell}^{m}(u)=\left(\begin{array}{cccc}{P*{\ell}^{m, 0}(u)} & {0} & {0} & {0} \ {0} & {P\_{\ell}^{m,+}(u)} & {P\_{\ell}^{m,-}(u)} & {0} \ {0} & {P\_{\ell}^{m,-}(u) P\_{\ell}^{m,+}(u)} & {0} \ {0} & {0} & {0} & {P\_{\ell}^{m, 0}(u)}\end{array}\right) $$ with $$P\_{\ell}^{m, \pm}(u)=\frac{1}{2}\left\[P\_{\ell}^{m,-2}(u) \pm P\_{\ell}^{m, 2}(u)\right]$$. The functions $$P\_{\ell}^{m, 0}(u)$$ and $$P\_{\ell}^{m, \pm 2}(u)$$ are the generalized spherical functions. We define $$P\_{\ell}^{m,n}(u)=0 \quad \text { for } \ell<\max {|m|,|n|}$$. For $$\ell>\max {|m|,|n|}$$ the functions are calculated as follow ($$i$$ is the imaginary unit). $$ \begin{aligned} P\_{m}^{m, 0}(u)=(2 i)^{-m}\left\[\frac{(2 m) !}{m ! m !}\right]^{1 / 2}\left(1-u^{2}\right)^{m / 2} \quad(m \geq 0) \end{aligned} \\ \begin{aligned} P\_{2}^{0, \pm 2}(u)=\frac{-1}{4} \sqrt{6}\left(1-u^{2}\right) \quad(m=0)\end{aligned} \\ \begin{aligned} P\_{2}^{1, \pm 2}(u)=\pm(2 i)^{-1}\left(1-u^{2}\right)^{1 / 2}(1 \pm u) \quad(m=1)\end{aligned} \\ \begin{aligned} P\_{m}^{m, \pm 2}(u)=&-(2 i)^{-m}\left\[\frac{(2 m) !}{(m+2) !(m-2) !}\right]^{1 / 2} \times(1-u)^{(m \mp 2) / 2}(1+u)^{(m \pm 2) / 2} \quad(m \geq 2) \end{aligned} $$ and the recurrence relations $$ \sqrt{(\ell+1)^{2}-m^{2}} P\_{\ell+1}^{m, 0}(u)=(2 \ell+1) u P\_{\ell}^{m, 0}(u)-\sqrt{\ell^{2}-m^{2}} P\_{\ell-1}^{m, 0}(u) $$ The equations above are taken from de Hann's 1987 paper, `eq. 74 to 81`, in $$\textit{The adding method for multiple scattering calculations of polarized light}$$. [http://articles.adsabs.harvard.edu//full/1987A%26A...183..371D/0000378.000.html](http://articles.adsabs.harvard.edu/full/1987A%26A...183..371D/0000378.000.html) The rotations of reference plane are implicitly accounted for in the expansion method. In the following code, `line 3~line6` uses functions defined outside. `plm0` calculates cases with m=0 and n=0, `plm2` calculates cases with m=2 and n=2, `plmn2` calculates cases with m=2, n=-2, and `plm2` calculates cases with m=0 and n=2. ```python nterms = length(data2); p00 = plm0(0,nterms-1,ctheta); % m=0, n=0 p22 = plm2(2,nterms-1,ctheta); % m=2, n=2 p2n2 = plmn2(2,nterms-1,ctheta); %P subscript m=2, n=-2 p02 = plm2(0,nterms-1,ctheta); %P subscript m=0, n=2 ``` Read in the phase moments from the `'mie_output_small2.dat'` file to alphas and betas (`line 2~8`). Suppressing the $$a = b$$ dependence, the following equations are for cases of non-spherical particles. $$ \begin{aligned} a\_{1}(\Theta)=\sum\_{\ell=0}^{2 N-1} \alpha\_{1}^{\ell} P\_{\ell}^{0,0}(\cos \Theta) \\ a\_{2}(\Theta)+a\_{3}(\Theta)=\sum\_{\ell=2}^{2 N-1}\left(\alpha\_{2}^{\ell}+\alpha\_{3}^{\ell}\right) P\_{\ell}^{2,2}(\cos \Theta) \\ a\_{2}(\Theta)-a\_{3}(\Theta)=\sum\_{\ell=2}^{2 N-1}\left(\alpha\_{2}^{\ell}-\alpha\_{3}^{\ell}\right) P\_{\ell}^{2,-2}(\cos \Theta) \\ a\_{4}(\Theta) &=\sum\_{\ell=0}^{2 N-1} \alpha\_{4}^{l} P\_{\ell}^{0,0}(\cos \Theta) \\ b\_{1}(\Theta) &=\sum\_{\ell=2}^{2 N-1} \beta\_{1}^{\ell} P\_{\ell}^{0,2}(\cos \Theta) \\ b\_{2}(\Theta) &=\sum\_{\ell=2}^{2 N-1} \beta\_{2}^{\ell} P\_{\ell}^{0,2}(\cos \Theta) \end{aligned} $$ The general form of Stokes scattering matrix Fs (which is the ensemble-averaged Mueller matrix averaging over a small volume containing an ensemble of particles) is of the form: $$ \mathbf{F}*{S}(\Theta)=\left(\begin{array}{cccc}{a*{1}(\Theta)} & {b\_{1}(\Theta)} & {0} & {0} \ {b\_{1}(\Theta)} & {a\_{2}(\Theta)} & {0} & {0} \ {0} & {0} & {a\_{3}(\Theta)} & {b\_{2}(\Theta)} \ {0} & {0} & {-b\_{2}(\Theta)} & {a\_{4}(\Theta)}\end{array}\right) $$ Note that there are six independent elements for non-spherical particles. ```python for i = 1:1 %1:numwl data2(:, :, i) = rawdata2((1+nmom*(i-1):nmom*i), 1:end); alpha1(:, i, 2) = data2(:, 1, i); alpha2(:, i, 2) = data2(:, 2, i); alpha3(:, i, 2) = data2(:, 3, i); alpha4(:, i, 2) = data2(:, 4, i); beta1(:, i, 2) = data2(:, 5, i); beta2(:, i, 2) = data2(:, 6, i); a1(:, i, 1) = alpha1(:,i, 2)' * p00; % Eq.68 a4(:, i, 1) = alpha4(:,i, 2)' * p00; % Eq.71 b1(:, i, 1) = beta1(:,i, 2)' * p02; % Eq.72 b2(:, i, 1) = beta2(:,i, 2)' * p02; % Eq.73 a2(:, i, 1) = .5*((alpha2(:,i, 2) + alpha3(:,i, 2))' * p22 + (alpha2(:,i, 2) - alpha3(:,i, 2))' * p2n2); %Eq.69+70 a3(:, i, 1) = .5*((alpha2(:,i, 2) + alpha3(:,i, 2))' * p22 - (alpha2(:,i, 2) - alpha3(:,i, 2))' * p2n2); %Eq.69+70 end ``` If we want to study spherical particles, there are four independent components, $$a\_{1}(\Theta), a\_{3}(\Theta), b\_{1}(\Theta), b\_{2}(\Theta)$$ , given $$a\_{2}(\Theta)=a\_{1}(\Theta), a\_{4}(\Theta)=a\_{1}(\Theta).$$ The code should be adjusted as follows $$ \mathbf{F}*{S}(\Theta)=\left(\begin{array}{cccc}{a*{1}(\Theta)} & {b\_{1}(\Theta)} & {0} & {0} \ {b\_{1}(\Theta)} & {a\_{1}(\Theta)} & {0} & {0} \ {0} & {0} & {a\_{3}(\Theta)} & {b\_{2}(\Theta)} \ {0} & {0} & {-b\_{2}(\Theta)} & {a\_{3}(\Theta)}\end{array}\right) $$ ```python for i = 1:1 %1:numwl data2(:, :, i) = rawdata2((1+nmom*(i-1):nmom*i), 1:end); alpha1(:, i, 2) = data2(:, 1, i); alpha2(:, i, 2) = data2(:, 2, i); alpha3(:, i, 2) = data2(:, 3, i); alpha4(:, i, 2) = data2(:, 4, i); beta1(:, i, 2) = data2(:, 5, i); beta2(:, i, 2) = data2(:, 6, i); a1(:, i, 1) = alpha1(:,i, 2)' * p00; % Eq.68 a4(:, i, 1) = alpha4(:,i, 2)' * p00; % Eq.71 b1(:, i, 1) = beta1(:,i, 2)' * p02; % Eq.72 b2(:, i, 1) = beta2(:,i, 2)' * p02; % Eq.73 a2(:, i, 1) = a1(:,i, 1) a3(:, i, 1) = a4(:,i, 1) end ``` {% hint style="info" %} Note: In here, a more precise way to do it is to output the a1 to a4 and b1, b2 directly from the mie code. The method above, to compute them from alpha1\~alpha4 and beta1 beta2 is actually taking a detour. {% endhint %} So far, we have calculated the scattering matrix $$\mathbf{F}*{S}$$. The phase matrix $$\mathbf{P}*{S}$$ is related to the scattering matrix by $$\mathbf{P}*{S}\left(u^{\prime}, \phi^{\prime}, u, \phi\right)=\mathbf{R}*{S}\left(\pi-\sigma\_{2}\right) \mathbf{F}*{S}(\Theta) \mathbf{R}*{S}\left(-\sigma\_{1}\right)$$. $$\mathbf{R}\_{S}$$ is a rotation matrix used to rotate reference planes, and it transforms phase matrix from the scattering plane to local meridian plane of reference. $$ \mathbf{R}\_{S}(\eta)=\left\[\begin{array}{cccc}{1} & {0} & {0} & {0} \ {0} & {\cos (2 \eta)} & {-\sin (2 \eta)} & {0} \ {0} & {\sin (2 \eta)} & {\cos (2 \eta)} & {0} \ {0} & {0} & {0} & {1}\end{array}\right] $$ Then use $$\mathbf{P}*{C}=\mathbf{D}^{-1} \mathbf{P}*{S} \mathbf{D}$$ (`line 10` in the following code snippet). $$\mathbf{F}*{C}(\Theta)=\mathbf{D}^{-1} \mathbf{F}*{S}(\Theta) \mathbf{D}$$, where $$\mathbf{F}\_{C}(\Theta)$$ is the Chandrasekhar's Stocks vector representation. Another way to find $$\mathbf{P}*{C}$$ is by first finding the transformed rotation matrix elements $$\tilde{\mathbf{R}}*{S}\left(\pi-\sigma\_{2}\right)$$ and $$\tilde{\mathbf{R}}*{S}\left(-\sigma*{1}\right)$$, and $$\mathbf{P}*{C}=\tilde{\mathbf{R}}*{S} (\pi-\sigma\_{2})\mathbf{P}*{S} \tilde{\mathbf{R}}*{S}(-\sigma\_{1})$$**.** The radiative transfer equitation for diffuse polarized radiation, described in terms of the Stokes equations is given by $$ \mu \frac{d \vec{I}(\tau, \mu, \phi)}{d \tau}=\vec{I}(\tau, \mu, \phi)-\frac{a(\tau)}{4 \pi} \int\_{0}^{2 \pi} d \phi^{\prime} \int\_{-1}^{1} d \mu^{\prime} \vec{P}\left(\tau, \mu, \phi ; \mu^{\prime}, \phi^{\prime}\right) \vec{I}\left(\tau, \mu^{\prime}, \phi^{\prime}\right)+\vec{S}(\tau, \mu, \phi) $$ Ignore the multiple-scattering term, we have $$\mu \frac{d \vec{I}(\tau, \mu, \phi)}{d \tau}=\vec{I}(\tau, \mu, \phi)+\vec{S}(\tau, \mu, \phi)$$ The source term is defined as $$\vec{S}(\tau, \mu, \phi)=\frac{a(\tau) F^{s}}{4 \pi} \vec{P}\left(\mu^{\prime}, \phi^{\prime} ; \mu, \phi\right) e^{-\tau / \mu\_{0}}$$ The diffused intensity in upward direction can be derived as $$\vec{I}*{d}^{+}(\tau, \mu, \phi)=\frac{a(\tau) \mu*{0} F^{s}}{4 \pi\left(\mu+\mu\_{0}\right)} \vec{P}\left(\mu^{\prime}, \phi^{\prime} ; \mu, \phi\right)\left\[e^{-\tau / \mu\_{0}}-e^{-\left\[\left(\tau^{*}-\tau\right) / \mu+\tau^{*} / \mu\_{0}\right]}\right]$$ For the intensity at the top of the layer ( $$\tau=0$$ ), the above equation becomes $$ \vec{I}*{d}^{+}(0, \mu, \phi)=\frac{a(\tau) \mu*{0} F^{s}}{4 \pi\left(\mu+\mu\_{0}\right)} \vec{P}\left(\mu^{\prime}, \phi^{\prime} ; \mu, \phi\right)\left\[1-e^{-\left(\tau^{*} / \mu+\tau^{*} / \mu\_{0}\right)}\right] $$ Isolate the Azimuthal dependence, we get the Chandraskhar's representation of $$\vec{I}\_{d}^{+}(\tau, \mu, \phi)$$ $$ \begin{aligned} \mathbf{I}(\tau, u, \phi)=& \sum\_{m=0}^{2 N-1}\left{\mathbf{I}*{c}^{m}(\tau, u) \cos m\left(\phi*{0}-\phi\right)+\mathbf{I}*{s}^{m}(\tau, u) \sin m\left(\phi*{0}-\phi\right)\right} \ \mathbf{Q}*{\delta}(\tau, u, \phi)=& \sum*{m=0}^{2 N-1}\left{\mathbf{Q}*{c \delta}^{m}(\tau, u) \cos m\left(\phi*{0}-\phi\right)+\mathbf{Q}*{s \delta}^{m}(\tau, u) \sin m\left(\phi*{0}-\phi\right)\right} \end{aligned} $$ As for the Stocks representation, we just need to read out the $$\vec{I}\_{d}^{+}(0, \mu, \phi)$$ calculated using the equation from the above. $$ \mathbf{I}*{C}=\left\[I*{\ell}, I\_{r}, U, V\right]^{T}=\mathbf{I}*{S}=\left\[I*{ |}, I\_{\perp}, U, V\right]^{T} $$ ````python for j = 1:length(phi) for i = 1:length(theta) Fs(:,:,(j-1)*181+i) = [a1(i,1), b1(i,1), 0, 0; b1(i,1), a2(i,1), 0, 0; 0, 0, a3(i,1), b2(i,1); 0, 0, -b2(i,1), a4(i,1)]; [Rs1(:,:,(j-1)*181+i), Rs2(:,:,(j-1)*181+i)] = Rs(mup, ctheta(i), cTheta(i)); Fc(:,:,(j-1)*181+i) = D^-1 * Fs(:,:,(j-1)*181+i) * D; Ps(:,:,(j-1)*181+i) = Rs2(:,:,(j-1)*181+i) * Fs(:,:,(j-1)*181+i) * Rs1(:,:,(j-1)*181+i); Pc(:,:,(j-1)*181+i) = D^-1 * Ps(:,:,(j-1)*181+i) * D; Rs1_1(:,:,(j-1)*181+i) = D^-1 * Rs1(:,:,(j-1)*181+i) * D; Rs2_1(:,:,(j-1)*181+i) = D^-1 * Rs2(:,:,(j-1)*181+i) * D; Pc_1(:,:,(j-1)*181+i) = Rs2_1(:,:,(j-1)*181+i) * Fc(:,:,(j-1)*181+i) * Rs1_1(:,:,(j-1)*181+i); Id(:,i,j) = a*mup/4/pi/(abs(ctheta(i))+mup)*(1-exp(-(tau/abs(ctheta(i))+tau/mup)))*Pc(:,:,(j-1)*181+i)*F; Im(j,i) = 1*(real(Id(1,i,j))+real(Id(2,i,j))); Qm(j,i) = 1*(real(Id(1,i,j))-real(Id(2,i,j))); Il(j,i) = real(Id(1,i,j)); Ir(j,i) = real(Id(2,i,j)); U(j,i) = real(Id(3,i,j)); V(j,i) = real(Id(4,i,j)); end ``` some codes to plot the graphs ``` end ```` Now it's time to scale it to two-slab and multi-slab cases. ### Single Multi-layer Slab For a multilayered slab the corresponding solution is $$ \begin{aligned} \mathbf{I}*{p}^{+}(\tau, \mu, \phi)=& \frac{\mu*{0}}{\left(\mu\_{0}+\mu\right)} \sum\_{n=p}^{L}\left{\mathbf{X}*{n}^{+}(\mu, \phi)\right.\ &\left. \times\left\[e^{-\left\[\tau*{n-1} / \mu\_{0}+\left(\tau\_{n-1}-\tau\right) / \mu\right]}-e^{-\left\[\tau\_{n} / \mu\_{0}+\left(\tau\_{n}-\tau\right) / \mu\right]}\right]\right} \end{aligned} $$ with $$\tau\_{n-1}$$ replaced by $$\tau$$ for n=p. p is the layer at which level that we are trying to evaluate, $$\tau\_{p}$$is the optical depth from the top of the slab to the p-th layer, and n is the total number of layer. ### Single, Two-layer slab $$\tau$$ is the depth measured from the top slab to any layer p. The subscript numbers get larger as more layers are accumulated downwards. As an example, for a 2 layer slab, with $$\tau\_{1}=0.0025, \tau\_{2}=0.005, \text{and we are evaluating an optical depth at }\tau=0.0024, \text{which means }p=1, \text{the equations are}$$ $$ \begin{array}{l}{I\_{1}^{+}=\frac{\mu\_{0}}{\left(\mu\_{0}+\mu\right)} X\_{1}\left(e^{-\frac{\tau}{\mu{0}}}-e^{-\left(\frac{\tau\_{1}}{\mu\_{0}}+\frac{\tau\_{1}-\tau}{\mu}\right)}\right)} \ \frac{\mu\_{0}}{\left(\mu\_{0}+\mu\right)} {+X\_{2}\left(e^{-\left(\frac{\tau\_{1}}{\mu\_{0}}+\frac{\tau\_{1}-\tau}{\mu}\right)}-e^{-\left(\frac{\tau\_{2}}{\mu\_{0}}+\frac{\tau\_{2}-\tau}{\mu}\right)}\right)} \end{array} $$ If evaluate on the bottom layer, say $$\tau=0.0026, \text{which means }p=2$$ $$ I\_{2}^{+}=\frac{\mu\_{0}}{\left(\mu\_{0}+\mu\right)} X\_{2}\left(e^{-\frac{\tau}{\mu\_{0}}}-e^{-\left(\frac{\tau\_{2}}{\mu\_{0}}+\frac{\tau\_{2}-\tau}{\mu}\right)}\right) $$ Therefore, to achieve this functionality, just need to (1) Evaluate the range to know what p is\ (2) Loop from p to n ![2 Layer, tau=0.0026](https://492391472-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LiBoWdc5sh0EFrazIOe%2F-LpnqL2IEkCVwIOwAren%2F-Lpnrh0h6Q01TniN_JJP%2FIUV-multi.jpg?alt=media\&token=bcccb1c7-28a0-4325-ac0e-00e5236e623e) ![1 Layer, tau=0.0026](https://492391472-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LiBoWdc5sh0EFrazIOe%2F-LpnqL2IEkCVwIOwAren%2F-LpnrqCAezFbDomGET3o%2FIUV-single.jpg?alt=media\&token=b5ee7042-5eb0-4384-815e-7a3e4e017f6b) The results are the same. Similar results can be achieved