Constant calculated by Simpson's method is not giving correct result, whilst when i give direct numerical value, the code is runnung.

Hello everyone, I am stuck at this little error for a while now. Following is a code for population dynamics of a four level system, although I have almost completed the first step but I need to calculate INTEGRATED ABSORPTION CROSSECTION, from EFFECTIVE CROSSECTION, this is done by straightaway deviding the effe_cross by the value of integral in known limit (line 68-70). The value of integral is "0.76". Here in line 71,72,73 I am directly giving the numeric input manually but this works, and the same thing done in line 68,69,70 does not. Please suggest where i am wrong so that I don't have to put numeric input.

The correct output given an array of w12,w23,w34 as all non zero values, which can be seen when s12, s23, s34 are passed as direct constants and not calculated by comment part , I have tried to put the start and end part of this calculation.

It would be of great help. Thank you.

h1 = 6.034e-34;  %plank's constant
ti = 0; %initial time
tf = 40e-9;  %final time
h = 1e-12;  %step size which is taken very small as the differential equations are stiff, i have not given those here, and RK4 is not a good option 
n = abs((tf-ti)/h);  % no of points in time
time = linspace(ti,tf,n);



%speed of light m/s
c = 299792458;

%transition wavelength in meter
l12 = 540e-9;
l23 = 535e-9;
l34 = 618e-9;

%effective absorption crossections in cm^2 over 200MHz
% s_eff12 = 5.1e-13;
% s_eff23 = 3.5e-12;
% s_eff34 = 1.5e-14;
%calculation for integrated absorption crossection
fwhm_s12 = 0.2;    %reduced fwhm for first transition in ghz due to doppler broadened first level in sigma measurement
sigma_s12 = fwhm_s12/2.35;
fwhm_s23 = 0.2;    %fwhm for second transition due to doppler broadened first level in sigma measurement
sigma_s23 = fwhm_s23/2.35;
fwhm_s34 = 0.2;    %fwhm for third transition due to doppler broadened first level in sigma measurement
sigma_s34 = fwhm_s34/2.35;

% fu=@(x,sigma)((exp(-(((x)./sigma).^2)./2))./(sqrt(2.*pi).*sigma));  %here starts the calculation of integrated crossection
% a1=-0.1;
% b1=0.1;
% 
% m=100;
% h=(b1-a1)./m;
% 
% so_s12=0;
% se_s12=0;
% so_s23=0;
% se_s23=0;
% so_s34=0;
% se_s34=0;
% y_s12=zeros(1,m-1);
% y_s23=zeros(1,m-1);
% y_s34=zeros(1,m-1);
% x=zeros(1,m-1);
% for l=1:1:m-1
%     x(l)=a1+l*h;
%     y_s12(l)=fu(x(l),sigma_s12);
%      y_s23(l)=fu(x(l),sigma_s23);
%       y_s34(l)=fu(x(l),sigma_s34);
%     if rem(l,2)==1
%        so_s12=so_s12+y_s12(l);%sum of odd terms
%        so_s23=so_s23+y_s23(l);%sum of odd terms
%        so_s34=so_s34+y_s34(l);%sum of odd terms
%     else
%        se_s12=se_s12+y_s12(l); %sum of even terms
%        se_s23=se_s23+y_s23(l); %sum of even terms
%        se_s34=se_s34+y_s34(l); %sum of even terms
%     end
% end
% %final values of integrated crossection
% s112=(h/3*(gw(a1,sigma_s12)+gw(b1,sigma_s12)+4*so_s12+2*se_s12))
% s213=(h/3*(gw(a1,sigma_s23)+gw(b1,sigma_s23)+4*so_s23+2*se_s23))
% s314=(h/3*(gw(a1,sigma_s34)+gw(b1,sigma_s34)+4*so_s34+2*se_s34))
% s12=s_eff12./(h/3*(fu(a1,sigma_s12)+fu(b1,sigma_s12)+4*so_s12+2*se_s12))
% s23=s_eff23./(h/3*(fu(a1,sigma_s23)+fu(b1,sigma_s23)+4*so_s23+2*se_s23))
% s34=s_eff34./(h/3*(fu(a1,sigma_s34)+fu(b1,sigma_s34)+4*so_s34+2*se_s34))               %here ends the calculation for integrated crossection
s12=(5.1e-13)./0.76                   %direct numerical valued which by principle should not be given
s23=(3.5e-12)./0.76
s34=(1.5e-14)./0.76
%transition frequencies in hz
n12 = c/l12;
n23 = (c/l23);
n34 = (c/l34);

%power of lasers using area dairy date 10/03
p1 = 2*(27e-6);
p2 = 54e-6;
p3 = 700e-6;


n_sol = zeros(4,length(time)-1);
n_sol(:,1 )= [100,0,0,0];

gauss = @(t,m,sig,amp)amp.*(1./(sqrt(2*pi).*sig)).*exp(-((t-m).^2)./(2.*sig.^2));             %intesity gaussian in time
time(1) = 0;

%T=input('temperature in K T:');
T = 1900;
%M=input('mass in amu M:')
M = 176;


fwhm_i12 = 2;             %fwhm for intensity in frequency dye laser 1 in ghz
sigma_i12 = fwhm_i12/2.35;
fwhm_i23 = 2;             %fwhm for intensity in frequency dl2
sigma_i23 = fwhm_i23/2.35;
fwhm_i34 = 2;             %fwhm for intensity in frequency dl3
sigma_i34 = fwhm_i34/2.35;
% fwhm_s12 = 0.2;    %reduced fwhm for first transition in ghz due to doppler broadened first level in sigma measurement
% sigma_s12 = fwhm_s12/2.35;
% fwhm_s23 = 0.2;    %fwhm for second transition due to doppler broadened first level in sigma measurement
% sigma_s23 = fwhm_s23/2.35;
% fwhm_s34 = 0.2;    %fwhm for third transition due to doppler broadened first level in sigma measurement
% sigma_s34 = fwhm_s34/2.35;

gw_is=@(dlf,sigmal,sigmaa)((exp(-((dlf./sigmal).^2)./2)).*(exp(-((dlf./sigmaa).^2)./2)))./(2.*sigmal.*sigmaa.*pi);         




a = -4;
b = 4;
n_1=100;
% n_1 = input('Enter the  number of sub-intervals n: '); 
h_1 =(b-a)/n_1;
% if rem(n_1,2)==1
%     fprintf('\n Enter valid n!!!'); 
%     n=input('\n Enter n as even number ');
% end


i12 = zeros(1,length(time)-1);
i23 = zeros(1,length(time)-1);
i34 = zeros(1,length(time)-1);
integrated_is12 = zeros(1,length(time)-1);
integrated_is23 = zeros(1,length(time)-1);
integrated_is34 = zeros(1,length(time)-1);
y_is12 = zeros(1,length(time)-1);
y_is23 = zeros(1,length(time)-1);
y_is34 = zeros(1,length(time)-1);

area = 10e-2;    %in cm^2


for i=1:length(time)-1
    time(i+1)=time(i)+h;
    
    amp12 = p1./(20*h1*n12*area*6e-9);
    m12 = 20e-9;
    sig12 = 2.553e-9;
    amp23 = p2./(20*h1*n23*area*6e-9);
    m23 = 22e-9;
    sig23 = 2.553e-9; 
    amp34 = p3./(20*h1*n34*area*6e-9);
    m34 = 24e-9;
    sig34=2.553e-9;
    i12(i) = gauss(time(i),m12,sig12,amp12);
    i23(i) = gauss(time(i),m23,sig23,amp23);
    i34(i) = gauss(time(i),m34,sig34,amp34);
    
    
    so_is12 = 0;
    so_is23 = 0;
    so_is34 = 0;
    se_is12 = 0;
    se_is23 = 0;
    se_is34 = 0;
    for k = 1:n_1-1
        dlf(k) = a+k*h_1;
        y_is12(k) = gw_is(dlf(k),sigma_i12,sigma_s12);
        y_is23(k) = gw_is(dlf(k),sigma_i23,sigma_s23);
        y_is34(k) = gw_is(dlf(k),sigma_i34,sigma_s34);
        if rem(k,2)==1
            so_is12 = so_is12+y_is12(k);%sum of odd terms
            so_is23 = so_is23+y_is23(k);%sum of odd terms
            so_is34 = so_is34+y_is34(k);%sum of odd terms
        else
            se_is12 = se_is12+y_is12(k); %sum of even terms
            se_is23 = se_is23+y_is23(k); %sum of even terms
            se_is34 = se_is34+y_is34(k); %sum of even terms
        end
        integrated_is12(i) = (h_1/3*(gw_is(a,sigma_i12,sigma_s12)+gw_is(b,sigma_i12,sigma_s12)+4*so_is12+2*se_is12)).*(i12(i).*s12.*1e-9);
        integrated_is23(i) = (h_1/3*(gw_is(a,sigma_i23,sigma_s23)+gw_is(b,sigma_i23,sigma_s23)+4*so_is23+2*se_is23)).*(i23(i).*s23.*1e-9);
        integrated_is34(i) = (h_1/3*(gw_is(a,sigma_i34,sigma_s34)+gw_is(b,sigma_i34,sigma_s34)+4*so_is34+2*se_is34)).*(i34(i).*s34.*1e-9);
    end
end


w12 = (integrated_is12);
w23 = (integrated_is23);
w34 = (integrated_is34);

Answer 1

[Ishu]

Hi Vipul,

I tried to reproduce the error at my end and tried to run the code that you have provided.

In your code for the calculation of 'integrated crosssection' you are calculating 's112', 's213' and 's314' right above the 's12', 's23' and 's34'. And I have gone through your code, you have not used 's112', 's213' and 's314' anywhere in the code.

s112=(h/3*(gw(a1,sigma_s12)+gw(b1,sigma_s12)+4*so_s12+2*se_s12))  % you can  
s213=(h/3*(gw(a1,sigma_s23)+gw(b1,sigma_s23)+4*so_s23+2*se_s23))  % remove these
s314=(h/3*(gw(a1,sigma_s34)+gw(b1,sigma_s34)+4*so_s34+2*se_s34))  % 3 lines

s12=s_eff12./(h/3*(fu(a1,sigma_s12)+fu(b1,sigma_s12)+4*so_s12+2*se_s12))
s23=s_eff23./(h/3*(fu(a1,sigma_s23)+fu(b1,sigma_s23)+4*so_s23+2*se_s23))
s34=s_eff34./(h/3*(fu(a1,sigma_s34)+fu(b1,sigma_s34)+4*so_s34+2*se_s34))    

And to calculate these variables('s112', s'213', 's314') you are using 'gw()' function which is not defined in the code because of which you are getting the error. You can simply remove these lines of code and then try to run the code. After that I guess this will not produce any error and you will be able to calculate the results correctly with both the methods.

You have used function handles in your code, if you want to learn more about these handles you can refer these documentations:

https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html

https://www.mathworks.com/help/matlab/function-handles.html

Hope it helps!