Curve fitting with complexe custom function

I want to fit some experimental data to a custom equation of this type:

Where a(T) and c(T) are some complicated functions of T (T is my x-axis parameter).

I calculate a(T) from a for loop. In my code a(T) and c(T) are represented by preExpPu(i,j,k) and preExpU(i,j,k). So basically the equation from the Figure above is my equation represented in the for loop by condTot(i,j,k). I know the values of a(T) and c(T) because I have calculated them. However, what I would like to fit is the values of b and d.... (in the exponential functions). The difficult part is that a(T) and c(T) are functions of T so I could not fit the whole data to a function of a type (a/T) exp(-b/T) + (c/T) exp(-d/T) because a and b are simply not constants, by vary with T... Is there any way to fit directly in the for loop do get the values of b and d?

for i=1:numel(T)

for j=1:numel(x)

for k=1:numel(y)

K(i)=exp((-H/(k1*T(i)))+1);

%Calcul de la fréquence moyenne des phonons (s-1)

v(j,k)=((2*pi*k2)/h)*TempDeb(j,k);

%Calcul des concentrations en trous/electrons

p(i,j,k)=-x(j)+(K(i)-sqrt(4*K(i)*x(j)*(K(i)*(x(j)-2*y(k)+1)+2*(y(k)-x(j))-1)+4*K(i)*y(k)*((y(k)-1)*(K(i)-1))+(K(i)^2)+4*(x(j)^2)))/(2*(K(i)-1));

n(i,j,k)=p(i,j,k)+2*x(j);

% % %Calcul des conductivité électriques dues à Pu et à U (S/m)

preExpPu(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))));

preExpU(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))));

condPu(i,j,k)= (preExpPu(i,j,k))/(T(:,i))*exp(-c/(k1*(T(:,i))));

condU(i,j,k)= (preExpU(i,j,k))/(T(:,i))*exp(-d/(k1*(T(:,i))));

% % %Calcul de la conductivité électrique totale (S/m)

condTot(i,j,k)=condPu(i,j,k)+condU(i,j,k);

end

end

end

[Mathieu NOE]

hello

maybe you should share some data wit a working code so we can (try) help you better

[Plamen Bonev]

x=[0, 0.01, 0.03, 0.04, 0.05, 0.07];

% x=[0.08, 0.09, 0.1, 0.15, 0.16, 0.17];

% x=[1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2];

T=[500 : 50 : 3223]

y=[0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];

Na=6.02e23;

R=8.314;

%Enthalpie de réaction en eV

H=1.5;

%h en J.s

h=6.63e-34

%k en eV/K

k1=8.6173303e-5;

%k en J/K

k2=1.38e-23;

%Charge élémentaire en Coulombs

e=1.6e-19;

%Energies de migration en eV

EhPu=0.8;

EhU=0.3;

%rayons en Angstrom (m)

rU4=0.1001e-9;

rPu4=0.096e-9;

rPu3=0.11e-9;

rO2=0.1368e-9;

rOv=0.1367e-9;

%Masses atomiques en Amu

Mu=238;

Mpu=244;

Mo=15.9;

eps=25.85;

format shortE

p=zeros(length(T),length(x),length(y));

n=zeros(length(T),length(x),length(y));

K=zeros(length(T),1);

condPu=zeros(length(T),length(x),length(y));

condU=zeros(length(T),length(x),length(y));

condTot=zeros(length(T),length(x),length(y));

condElectron=zeros(length(T),length(x),length(y));

yU4=zeros(length(y),1);

yPu4=zeros(length(x), length(y));

yPu4stoech=zeros(length(y),1);

yPu3=zeros(length(x),1);

yO2=zeros(length(x),1);

yO2stoech=(2/3);

yOv=zeros(length(x),1);

rmeanCat=zeros(length(x), length(y));

rmeanAn=zeros(length(x),1);

a=zeros(length(x),length(y));

vl=zeros(length(x),length(y));

vt=zeros(length(x),length(y));

N=zeros(length(x),1);

TempDeb=zeros(length(x),length(y));

v=zeros(length(x),length(y));

VaMetre=zeros(length(x),length(y));

VaAngstr=zeros(length(x),length(y));

Mmoy=zeros(length(y),1);

rmoy=zeros(length(y),1);

SommeGamma=zeros(length(x),length(y));

V=zeros(length(x),length(y));

vp=zeros(length(x),length(y));

C=zeros(length(x),length(y));

A=zeros(length(x),length(y));

LTE=zeros(length(x),length(y));

a0=zeros(length(x),length(y));

a1=zeros(length(x),length(y));

a2=zeros(length(x),length(y));

a3=zeros(length(x),length(y));

alpha=zeros(length(x),length(y));

MmoyMOXkG=zeros(length(x),length(y));

MmoyMOXamu=zeros(length(x),length(y));

rhoTh=zeros(length(x),length(y));

BulkModulus=zeros(length(x),length(y));

Cv=zeros(length(x),1);

Vm=zeros(length(x),length(y));

Grun=zeros(length(x),length(y));

constB=zeros(length(x),length(y));

B=zeros(length(x),length(y));

CondPhonique=zeros(length(T),length(x),length(y));

preExpPu=zeros(length(T),length(x),length(y));

preExpU=zeros(length(T),length(x),length(y));

% %y=0.2

for i=1:numel(T)

for j=1:numel(x)

for k=1:numel(y)

%Pour le calcul du paramètre de maille

yU4(k)=(1-y(k));

yPu4(j,k)=(y(k)-2*x(j));

yPu4stoech(k)=y(k)/3;

yPu3(j)=(2*x(j));

yO2(j)=(2-x(j))/2;

yOv(j)=x(j)/2;

rmeanCat(j,k)=yPu3(j)*rPu3+yPu4(j,k)*rPu4+yU4(k)*rU4;

rmeanAn(j)=yO2(j)*rO2;

%Pour le calcul de Gamma(i)

rmoy(k)=(yU4(k)/3)*rU4+yPu4stoech(k)*rPu4+(yO2stoech)*rO2;

Mmoy(k)=(yU4(k)/3)*Mu+yPu4stoech(k)*Mpu+(yO2stoech)*Mo;

%Pour le calcul de la fréquence moyenne des phonons

vl(j,k)=5358*(1-0.7279*x(j))*(1+0.040*y(k));

vt(j,k)=2750*(1-1.0545*x(j))*(1+0.043*y(k));

N(j)=4*(3-x(j));

%Calcul du paramètre de maille

a(j,k)=(4/(3^(1/2)))*(rmeanCat(j,k)+rmeanAn(j));

% Calcul du volume atomique moyen

VaMetre(j,k)=((a(j,k))^3)/N(j); %en m^3

VaAngstr(j,k)=((a(j,k)*10^(10))^3)/N(j); %en A^3

%Calcul du volume atomique moyen par molecule de MOX

V(j,k)=((a(j,k))^3)/4;

% Calcul de la somme de Gamma(i);

SommeGamma(j,k)=(((x(j))/3)*(((-Mo^2)/(Mmoy(k))^2)+(eps/(rmoy(k)^2))*(2*(rPu3^2-rPu4^2)-(rOv^2-rO2^2))))+(0.035/0.886);

%

%Calcul de la temperature de Debye (K)

TempDeb(j,k)=(h/k2)*((9*N(j))/(4*pi*(a(j,k)^3)*(1/(vl(j,k))^3+2/(vt(j,k))^3)))^(1/3);

%Calcul de la vitesse moyenne des phonons (m/s)

vp(j,k)=((2*pi*k2*TempDeb(j,k))/(h))*((V(j,k))/(6*(pi^2)))^(1/3);

%Calcul de C

C(j,k)=((pi^2)*VaMetre(j,k)*TempDeb(j,k))/(3*(vp(j,k)^2)*h);

%Calcul de A

A(j,k)=C(j,k)*SommeGamma(j,k);

%--------------------------------------------------------------------------

%Pour le calcul du LTE

a0(j,k)=(-2.8809e-3)+((0.0301e-3)*y(k))-(4.3954e-3)*x(j)+(0.0156e-3*((y(k))^2))-(15.1759e-3)*((x(j))^2)+(2.5642e-3)*y(k)*x(j);

a1(j,k)=(9.5024e-6)+((-0.1864e-6)*y(k))+((15.8173e-6)*x(j))+((-0.0229e-6)*(y(k))^2)+((7.6258e-6)*(x(j))^2)+((-7.5789e-6)*y(k)*x(j));

a2(j,k)=(2.0894e-10)+(2.9483e-10)*y(k)+(-19.9227e-10)*x(j)+(-1.0355e-10)*(y(k))^2+(73.8931e-10)*(x(j))^2+(11.6442e-10)*y(k)*x(j);

a3(j,k)=(4.4096e-13)+(-1.4263e-13)*y(k)+(23.5638e-13)*x(j)+(0.0251e-13)*(y(k))^2+(-54.751e-13)*(x(j))^2+(-14.418e-13)*y(k)*x(j);

%Calcul de LTE à T=800K (pour B)

LTE(j,k)=a0(j,k)+a1(j,k)*800+a2(j,k)*(800^2)+a3(j,k)*(800^3);

%Calcul du coefficient volumétrique d'expansion thermique alpha (pour B)

alpha(j,k)=3*((LTE(j,k))/(800-300));

%Calcul de la densité théorique (pour B)

MmoyMOXkG(j,k)=((1-y(k))/3)*(Mu*10^(-3))+(y(k)/3)*(Mpu*10^(-3))+((2-x(j))/3)*(Mo*10^(-3)); %en kg/mol

MmoyMOXamu(j,k)=((1-y(k))/3)*(Mu)+(y(k)/3)*(Mpu)+((2-x(j))/3)*(Mo); %en g/mol

rhoTh(j,k)=(4*(3-x(j))*(MmoyMOXkG(j,k)))/((a(j,k)^3)*Na);

%Calcul du Bulk Modulus (Pa) (pour B)

BulkModulus(j,k)=((rhoTh(j,k))/3)*(3*(vl(j,k)^2)-4*(vt(j,k))^2);

%Calcul du Cv (pour B)

Cv(j)=3*(3-x(j))*R;

%Calcul du volume molaire (pour B)

Vm(j,k)=((a(j,k)^3)*Na)/(4);

%Calcul du paramètre de Gruneisen (pour B)

Grun(j,k)=(alpha(j,k)*Vm(j,k)*BulkModulus(j,k))/(Cv(j));

%Calcul de la constante B' (pour B)

constB(j,k)=(2.43e-8)/(1-(0.514/Grun(j,k))+(0.228/Grun(j,k)^2));

%Calcul de B

B(j,k)=((1/constB(j,k))*((Grun(j,k)^2)*(4*(3-x(j)))^(2/3))/(MmoyMOXamu(j,k)*(TempDeb(j,k)^3)*VaAngstr(j,k)^(1/3)))*(10^(-2));

%Calcul de la conductivité phonique

CondPhonique(i,j,k)=1/(A(j,k)+B(j,k)*T(i));

%---------------------------------------------------------------------------

%Calcul de la constante d'équlibre de la réaction de disproportiatiation

K(i)=exp((-H/(k1*T(i)))+1);

%Calcul de la fréquence moyenne des phonons (s-1)

v(j,k)=((2*pi*k2)/h)*TempDeb(j,k);

%Calcul des concentrations en trous/electrons

p(i,j,k)=-x(j)+(K(i)-sqrt(4*K(i)*x(j)*(K(i)*(x(j)-2*y(k)+1)+2*(y(k)-x(j))-1)+4*K(i)*y(k)*((y(k)-1)*(K(i)-1))+(K(i)^2)+4*(x(j)^2)))/(2*(K(i)-1));

n(i,j,k)=p(i,j,k)+2*x(j);

% % %Calcul des conductivité électriques dues à Pu et à U (S/m)

condPu(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))))*exp(-EhPu/(k1*(T(:,i))));

condU(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))))*exp(-EhU/(k1*(T(:,i))));

% preExpPu(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))));

% preExpU(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))));

% condPu(i,j,k)=((preExpPu(i,j,k))/(T(:,i)))*exp(-EhPu/(k1*(T(:,i))));

% condU(i,j,k)=((preExpU(i,j,k))/(T(:,i)))*exp(-EhU/(k1*(T(:,i))));

% % %Calcul de la conductivité électrique totale (S/m)

condTot(i,j,k)=condPu(i,j,k)+condU(i,j,k);

% condTot(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*(n(i,j,k)+p(i,j,k))*(1-n(i,j,k)-p(i,j,k)))*(exp(-0.1/(k1*(T(:,i)))));

% Preexp=((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*(n(i,j,k)+p(i,j,k))*(1-n(i,j,k)-p(i,j,k))

%Calcul de la conductivité électronique (W/m/K)

condElectron(i,j,k)=((k2/e)^2)*(T(:,i))*(condPu(i,j)*condU(i,j)/condTot(i,j))*((H/(k1*T(:,i)))^2);

end

end

end