The file 'thermo_refractive_COMSOL_1D_axisymmetric.m' found in the SVN repository https://nodus.ligo.caltech.edu:30889/svn/trunk/comsol/thermo-refractive/ performs the data extraction from the COMSOL simulation
and computes quantities as given in Heinert's paper.
This is the where the dissipation is calculated from the temperature gradient and the use of FDT is made to evaluate the linearized PSD.
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
The first few lines of the code extracts the variable values out of the structure 'param' and stores it in new variables for the ease of coding.
The important portion comes in with the iterative structure (for loop) in line 29.
A point wise summary of the steps done would be as follows -
* The finite elements in space and time are defined in the arrays r0 and t0 according to the time and radial slice steps dt and dR until R, the radius and the t_end, the time until which simulation runs
* A matrix called dat is created which is used to store the values of the temperature gradient as returned by COMSOL.
It stores the values as a function of r along the columns and as a function of t along the rows.
Thus, moving vertically down the matrix in a column would give the values of the temperature gradient for a fixed r as a function of t
* Two other row vectors datrabs and datrphs having lengths of the array r0 are used to store the values of the square of the Fourier coefficients and the phase angle.
* The extraction of data from COMSOL is done in the lines 71-74
for i2=1:length(t0)
tmp = mphinterp(model,'ht.gradTr','coord',r0,'T',t0(i2));
dat(:,i2) = tmp;
end
* After the above step, the matrix dat is filled with the values of the temperature gradient as mentioned previously.
* Next to evaluate to the Fourier coefficients, the second half of the simulation time is used. Probably to let the transients die down as mentioned in the ELOG post by the author.
* To find out the Fourier components, we multiply the desired function by sin(omega*t) or cos(omega*t) and integrate over the period of the function. There is also a prefactor of 1/L where L is the length of the period
over which the function is defined.
* To do the above exercise, the sine and cosine of omega*t is separately evaluated and stored in skernel and ckernel. Note that the time interval as mentioned above is the second half of the simulation time.
* Now for all the radial slices, the Fourier coefficients are extracted using the statements in for loop in line 83
for i2=1:length(r0)
tseries=detrend(dat(i2,(length(t0)+1)/2:length(t0)));
tmp1=mean(tseries.*skernel);
tmp2=mean(tseries.*ckernel);
datrabs(i2)=2*sqrt(tmp1^2+tmp2^2);
datrphs(i2)=atan2(tmp2,tmp1)/pi*180;
end
Using the definition of mean of a function ' f ' in the continuous case as f_mean = integral{ f(t)dt} / integral{ dt } , over the values of t
tmp1 has used ' f ' as tseries.*skernel which means that one multiplies the time signal with sin(omega*t) and integrates. tmp2 has used cosine instead
Thus, tmp1 and tmp2 are the Fourier coefficients of tseries at the frequency fmod defined right close to the first for loop in line 29.
* The coefficients are squared and added in the line 88
datrabs(i2)=2*sqrt(tmp1^2+tmp2^2);
This quantity how much of the frequency ' fmod' is present in the time signal.
datrphs stores the same for the phase shift phi. Note however, all of this is done as a function of r. The plotting is done in the following lines. Note that while running the code, as a result, the plots are seen to
get updated each time the loop [line 29] iterates.
* Following this, is the part where W_diss, the dissipation is calculated using formula [The expression used does not match Eq. 15 in Heinert's. Any comments regarding the formula used for W_diss are welcome].
The heinert's eq(15) goes as - W_diss = pi*H*kappa/T0 * integral{ grad_T^2 * rdr }
The quantity datrabs is squared and the integral mentioned above is calculated w.r.t. r using the trapz algorithm from Matlab.
The value for each frequency i.e. each 'fmod' is converted to a string and displayed on screen.
* The reason why the grad_T was changed to its frequency domain before integration is suspected to be something related to the Parseval's Theorem. However, more details or correction on the same is welcome. |