Improved the function cpNb3Sn_Arp_Knapp_mat.m: 1) it now works when Tc, Tc0,...

Improved the function cpNb3Sn_Arp_Knapp_mat.m: 1) it now works when Tc, Tc0, Tcs are either scalars or vectors. 2) if Tc, Tc0, and Tcs are all vectors with identical elements. If so, they'll be changed to scalars to speed up calculation. 3) Small speed improvement by commenting redundant "ps(idxAboveTc) = 0;". 4) changed code reformatting to make it more readable. NOTE: The edited function is expected to return the same values as before (when Tc, Tc0, Tcs are scalars).

models/matlab/cpNb3Sn_Arp_Knapp_mat.m

@@ -6,6 +6,15 @@ function cpNb3Sn = cpNb3Sn_Arp_Knapp_mat(T,Tc,Tc0,Tcs)
 % C #            superconductors obtained from heat capacity measurements,
 % C #            Phys.Rev B, 13, no.9, pp 3783-3789, 1976.
 
+
+% Check whether Tc, Tc0, and Tcs are all vectors with identical elements
+% If so, they'll be changed to scalars to speed up calculation
+if all(Tc == Tc(1)) && all(Tc0 == Tc0(1)) && all(Tcs == Tcs(1))
+    Tc = Tc(1);
+    Tc0 = Tc0(1);
+    Tcs = Tcs(1);
+end
+
 pn=zeros(size(T));
 ps=zeros(size(T));
 p=zeros(size(T));
@@ -13,10 +22,10 @@ p=zeros(size(T));
 density = 8950.0;
 
 % Normal component
-idx1=find(T<=10);
-idx2=find(T>10&T<=20);
-idx3=find(T>20&T<=400);
-idx4=find(T>400);
+idx1 = find(T <= 10);
+idx2 = find(T > 10 & T <= 20);
+idx3 = find(T > 20 & T <= 400);
+idx4 = find(T > 400);
 
 AA  =   38.2226876;
 BB  =   -848.36422;
@@ -31,35 +40,49 @@ nb  =   2;
 nc  =   3;
 nd  =   4;
 
-pn(idx1)=(7.5475e-3)*T(idx1).^2;
-pn(idx2)=(-0.3+0.00375*T(idx2).^2)/0.09937;
-pn(idx3)=AA*T(idx3)./(a+T(idx3)).^na+BB*T(idx3).^2./(b+T(idx3)).^nb + CC*T(idx3).^3./(c+T(idx3)).^nc+DD*T(idx3).^4./(d+T(idx3)).^nd;
-pn(idx4)=250.8246;
+pn(idx1) = (7.5475e-3) * T(idx1).^2;
+pn(idx2) = (-0.3 + 0.00375 * T(idx2).^2) / 0.09937;
+pn(idx3) = AA * T(idx3) ./ (a + T(idx3)).^na + BB * T(idx3).^2 ./ (b + T(idx3)).^nb ...
+          + CC * T(idx3).^3 ./ (c + T(idx3)).^nc + DD * T(idx3).^4 ./ (d + T(idx3)).^nd;
+pn(idx4) = 250.8246;
 
 
 % Superconducting component
-idxBelowTc=find(T<Tc);
-idxBelowTcs=find(T<Tcs);
-idxAboveTc=find(T>=Tc);
+idxBelowTc  = find(T < Tc);
+idxBelowTcs = find(T < Tcs);
+idxAboveTc  = find(T >= Tc);
 
-TT=Tc/Tc0;
-dBc2Dt=-0.46306-0.067830*Tc;
-delcp=1500*(dBc2Dt^2)/(2*(27.2/(1+0.34*TT))^2-1);
-if Tc<=10
-    cpntc=7.5475e-3*Tc^2;
-elseif Tc<=20
-    cpntc=(-0.3+0.00375*Tc^2)/0.09937;
-end
+TT = Tc ./ Tc0;
+dBc2Dt = -0.46306 - 0.067830 * Tc;
+delcp = 1500 * (dBc2Dt.^2) ./ (2 * (27.2 ./ (1 + 0.34 * TT)).^2 - 1);
 
-ps(idxBelowTc)=(cpntc+delcp)*(T(idxBelowTc)/Tc).^3;
-ps(idxAboveTc)=0;
+if isscalar(Tc) % Calculation is the same, but code is slightly differnet if Tc is scalar or not
+    if Tc <= 10
+        cpntc = 7.5475e-3 * Tc^2;
+    elseif Tc <= 20
+        cpntc = (-0.3 + 0.00375 * Tc^2) / 0.09937;
+    else
+        cpntc = 0;
+    end
+    ps(idxBelowTc) = (cpntc + delcp) * (T(idxBelowTc) / Tc).^3;
+    F = (T(idxBelowTc) - Tcs) / (Tc - Tcs);
+else
+    cpntc = zeros(size(Tc));
+    idxTcLow = find(Tc <= 10);
+    idxTcMid = find(Tc > 10 & Tc <= 20);
+    cpntc(idxTcLow) = 7.5475e-3 * Tc(idxTcLow).^2;
+    cpntc(idxTcMid) = (-0.3 + 0.00375 * Tc(idxTcMid).^2) / 0.09937;
+    ps(idxBelowTc) = (cpntc(idxBelowTc) + delcp(idxBelowTc)) .* (T(idxBelowTc) ./ Tc(idxBelowTc)).^3;
+    F = (T(idxBelowTc) - Tcs(idxBelowTc)) ./ (Tc(idxBelowTc) - Tcs(idxBelowTc));
+end
+% ps(idxAboveTc) = 0; % No need for explicit assignment as ps is preallocated to zeros
 
-p(idxAboveTc)=pn(idxAboveTc);
-F=(T(idxBelowTc)-Tcs)/(Tc-Tcs);
-p(idxBelowTc)=F.*pn(idxBelowTc)+(1-F).*ps(idxBelowTc);
-p(idxBelowTcs)=ps(idxBelowTcs);
 
+% Calculate heat capacity
+p(idxAboveTc) = pn(idxAboveTc);
+p(idxBelowTc) = F .* pn(idxBelowTc) + (1 - F) .* ps(idxBelowTc);
+p(idxBelowTcs) = ps(idxBelowTcs);
 
-cpNb3Sn=density*p;
+cpNb3Sn = density * p;
 
-end
\ No newline at end of file
+end