How to sum up the functional sequence correctly?

I have a complex function that depends on 't' and 'n' F(n,t) I need to sum a series over 'n' and plot F(t).

Ef = 2.77*10^3;
Kb = physconst('boltzmann'); % 1.38*10^(-23)
double T;
%T = 0:366;
m = 9.1093837*10^(-31);
double Tc;
Tc = 1.2;
%T = 1.2; 
double t;
t = T/Tc;
int8 n;
%n = 0:1:366; %0:49;
double D;
%D = 10^(-8):10^(-5);
double ksi; % длина когерентности
%D = input('input D');
%ksi = input('input ksi');
D = 10^(-6);
ksi = 10^(-9);
d = D/ksi;
E = Ef/(pi*Kb*Tc);
hp = (626176*10^(-34))/(2*pi);
%complex p1;
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
p2 = sqrt(E+1i.*t.*(2.*n+1));
k0 = (ksi/hp)*sqrt(2.*m.*pi.*Kb.*Tc);
arg_of_lg = (1+exp(-1i*d*k0.*p1))/(1+exp(-1i*d*k0.*p2));
l_g = log(abs(arg_of_lg));
arg_1 = angle(1+exp(-1i*d*k0.*p1));
arg_2 = angle(1+exp(-1i*d*k0.*p2));
F = 0;
%F = 1/(2*n+1).*imag(l_g+1i*d*k0.*((p1-p2)./2)+1i*arg_1-1i*arg_2);
for n = 0:49
    %p1 = @Calculate_p1;
    %p2 = @Calculate_p2;
    F  = F + 1/(2*n+1).*imag(l_g+1i*d*k0.*((p1-p2)./2)+1i*arg_1-1i*arg_2);
end
F = -(1/d).*F;
F = F - (hp^2*pi^2*ksi)/(2*m*10^(-7));

[Jan]

What is the problem with the posted code?

Omit these useless lines:

double Tc;
double t;
int8 n;
double D;
double ksi;

You can simplify 9.1093837*10^(-31); to 9.1093837e-31. Note that the first is an expensive power operation while the second is a cheap constant.

Answer 1

[Torsten]

t = (100:1000).';
n = 0:49;
Ef = 2.77e3;
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837e-31;
Tc = 1.2;
D = 1e-6;
ksi = 1e-9;
d = D/ksi;
E = Ef/(pi*Kb*Tc);
hp = (626176*10^(-34))/(2*pi);
k0 = (ksi/hp)*sqrt(2.*m.*pi.*Kb.*Tc);
p1 = ((1+1i)/sqrt(2))*sqrt(t.*(2.*n+1));
p2 = sqrt(E+1i*t.*(2.*n+1));
arg_of_lg = (1+exp(-1i*d*k0.*p1))./(1+exp(-1i*d*k0.*p2));
l_g = log(abs(arg_of_lg));
arg_1 = angle(1+exp(-1i*d*k0.*p1));
arg_2 = angle(1+exp(-1i*d*k0.*p2));
F = 1./(2*n+1).*imag(l_g+1i*d*k0.*((p1-p2)./2)+1i*arg_1-1i*arg_2);
format long
F = sum(F,2)

F = 901×1

1.0e+10 * -1.047059455668692 -1.047059455668693 -1.047059455668693 -1.047059455668694 -1.047059455668695 -1.047059455668695 -1.047059455668696 -1.047059455668696 -1.047059455668696 -1.047059455668697

F = -(1/d).*F;
F = F - 1e7*(hp^2*pi^2*ksi)/(2*m);
plot(t,F)