% TE calculation using 2D axisymmetric model. This script calls the model
% file, extracts grad_T and does the calculations for the different
% frequency
% All parameter are defined in this script and passed to the COMSOL model
% via arguments
% Parameters, Quantitites in [m]
param.beam.radius = 0.05;
param.mirror.radius = 100 * param.beam.radius; 
%Define the radius of the mirror as 100 times that of the beam radius
param.mirror.height = 200 * param.beam.radius;
% Define the height of the cylinder to be 200 times the beam radius
param.COMSOL.total_radial_slices = 20;
param.COMSOL.total_z_slices = 40;

param.constant.boltzmann = 1.38065e-23;
param.constant.F0 = 1e-2;% Force amplitude
param.constant.T_amb = 300;
param.COMSOL.freq = 10.^[0:0.25:4];% frequency range
param.COMSOL.simulation_cycles = 10;% The number of cycles of Pressure Osc.
param.COMSOL.simulation_time_slice_per_cycle = 20; % Time steps in each cycle

% Parameters for fused silica
param.material = fused_silicadata;

% Assign values to different varibales
alpha = param.material.alpha;
E = param.material.E;
C = param.material.specific_heat;
kappa = param.material.thermal_conductivity;
sigma = param.material.sigma;
rho = param.material.density;
kB = param.constant.boltzmann;
T = param.constant.T_amb;
r0 = param.beam.radius;
F0 = param.constant.F0;
disp('Starting Time');
tic;
freq = param.COMSOL.freq;
psd = zeros(1,length(freq));% to store the PSD values

%%Calculation and analysis :

for i = 1:length(freq)
    param.COMSOL.fmod = freq(i);% COMSOL will be asked to simulate with this frequency
    dt = 1/(freq(i) * param.COMSOL.simulation_time_slice_per_cycle);
    param.COMSOL.simulation_time_step = dt;% time steps
    param.COMSOL.simulation_end_time = 1/freq(i) * param.COMSOL.simulation_cycles;% run simulation until this time
    t_end = param.COMSOL.simulation_end_time;
    disp(['Running COMSOL for frequency = ' num2str(freq(i))]);
    %% Call to COMSOL
    model = almost_infinite_2D_axisym_TE_4(param);
    %% Data Extraction
    % definitions
    r_slices = param.COMSOL.total_radial_slices;% assign new variables for convinience
    dR = param.mirror.radius/r_slices;% size of each R slice
    z_slices = param.COMSOL.total_z_slices;
    dz = param.mirror.height/z_slices;% size of each z slice
    t_slices = t_end/dt;
    disp('## Extracting data');
    % Store the data for each r, z for the last three cycles of simulation
    % time
    dat = zeros(r_slices,z_slices,floor(3/(freq(i)*dt)));% to store data from COMSOL for each r, z and 't'- for the last three cycles
    for m = 0:(r_slices-1)
        for k = 0:(z_slices-1)
            dat(m+1,k+1,:) = mphinterp(model,'ht.gradTmag','coord',[m*dR ; k*dz],'t', [(t_slices - floor(3/(freq(i)*dt)) + 1) : t_slices]);
            % The data is extracted for the last 3 cycles of simulation
        end
    end
    t_eval = dt * ((t_slices - floor(3/(freq(i)*dt)) + 1) : t_slices);% the time values for the last 3 cycles;
    temp = exp(j * 2*pi*freq(i) .* t_eval);% stores all exp(i * omega * t) values
    kernel = zeros(1,1,length(temp));
    kernel(1,1,:) = temp;% Now dimensions of kernel and dat_fc matches
    % Without the above modification, the assignment in the nested for
    % gives an error
    dat_fc = zeros(r_slices,z_slices);% matrix to store the fourier coefficients for each r, z time signal
    %% Fourier coefficients corresponding to fmod are extracted
    disp('## Extracting Fourier coefficients');
    for m = 1:r_slices
        for k = 1:z_slices
            dat_fc(m,k) = abs(2*mean(kernel.*dat(m,k,:)));% this is sqrt(an^2 + bn^2)
        end
    end
    disp('## Extraction Complete. Evaluating volume integral');
    %% Calculation of the integral of grad_T^2
    dat_fc = dat_fc.^2; % square the values since we need (an^2 + bn^2) instead of its sqrt
    dat_fc = dat_fc';% take a transpose so that going down a column you get all z values for a fixed r
    z = 0:dz:(param.mirror.height-dz);% the z values
    r = 0:dR:(param.mirror.radius-dR);% the r values
    z_integral = trapz(z,dat_fc);% integrate over the columns of the matrix dat_fc
    vol_integral = trapz(r,2*pi*r.*z_integral);% total volume integration of (an^2 + bn^2)
    % There must be an additional factor of 1/2 since the integral we want
    % is 1/2 (an^2 + bn^2)
    vol_integral = 0.5 * vol_integral;% Now we have time_avg <integral{ grad_T^2 dV}> as calculated in the paper
    psd(i) = vol_integral/(2*pi*freq(i))^2;% divide by omega^2 and store
    toc;
end


psd = 8*kB*kappa/F0^2 * psd; %Multiply all the prefactors in the formula of FDT
loglog(freq, psd,'b.');
hold on;
LT_constant = (8*(1 + sigma)^2*kappa*alpha^2*kB*T^2) / (sqrt(2*pi)*(rho*C)^2*r0^3);
loglog(freq, LT_constant./(2*pi.*freq).^2,'ro');
hold off;