Matrix close to singular or badly scaled problem in Gauss-Newton LMS calculations

I am trying to implement a function that minimizes the square of the residuals (LMS) for a nonlinear function. I am trying to use the Gauss-Newton method. The nonlinear function is the Fresnel reflection coeficients, which are calculated with another function I implemented, called bwTMM.

The minimizing function is called GNTMM. When I try the function with several different parameters i get the error:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In GNTMMs (line 58)

Any ideas? I insert 3 codes: The function bwTMM that calculates the Fresnel coefficients, the function that tries to fit the parameters, and the function in which I try the fitting function and get an error.

Thank you for taking the time to read:

function [GMs,GMp,phiout,Rs,Ts,Rp,Tp]=bwTMM(n,e,phi0,lambda)%e y lambda en mismas unidades

ko=2*pi/lambda;%calculo ko

phi=deg2rad(phi0);%defino el primer phi (angulo de incidencia)

Ms=eye(2);

Mp=eye(2);

for k=1:length(e)

phi=asin((n(k)/n(k+1))*sin(phi));

p=n(k+1)*cos(phi);%p es sqrt(epsilon/mu)cos(phi), pero si medio no magnetico mu=1 y, como sqrt(epsilon)=n, queda p

B=ko*n(k+1)*e(1)*cos(phi);

q=(1/n(k+1))*cos(phi);% Esto es para TM, o polarizacion paralela al plano de incidencia

%--------Matriz para polarizacion TE (plarizacion s)

ms11=cos(B);

ms12=-(1i/p)*sin(B);

ms21=-1i*p*sin(B);

ms22=cos(B);

Ms=Ms*[ms11,ms12;ms21,ms22];

%------Matriz polarizacion TM (p polarized)-------

mp11=cos(B);

mp12=-(1i/q)*sin(B);

mp21=-1i*q*sin(B);

mp22=cos(B);

Mp=Mp*[mp11,mp12;mp21,mp22];

end

phiout=rad2deg(asin((n(length(e)+1)/n(length(n)))*sin(phi)));

p1=n(1)*cos(phi0);

pl=n(end)*cos(phiout);

q1=(1/n(1))*cos(phi0);

ql=(1/n(end))*cos(phiout);

%-----Extraigo elementos de Ms global:-------

ms11g=Ms(1,1);

ms12g=Ms(1,2);

ms21g=Ms(2,1);

ms22g=Ms(2,2);

%------Extraigo elementos de Mp global:------

mp11g=Mp(1,1);

mp12g=Mp(1,2);

mp21g=Mp(2,1);

mp22g=Mp(2,2);

%----Coeficientes R Y T para polarizacion s------

%-----Calculo de r y t (coeficientes de amplitud)

dens=(ms11g+pl*ms12g)*p1+(ms21g+pl*ms22g);

rs=((ms11g+pl*ms12g)*p1-(ms21g+pl*ms22g))/dens;

ts=2*p1/dens;

%-----Coeficientes R y T (intensidad)-------

Rs=(abs(rs))^2;

Ts=(pl/p1)*(abs(ts))^2;

%----Coeficientes R Y T para polarizacion p------

%-----Calculo de r y t (coeficientes de amplitud)

denp=(mp11g+ql*mp12g)*q1+(mp21g+ql*mp22g);

rp=((mp11g+ql*mp12g)*q1-(mp21g+ql*mp22g))/denp;

tp=2*q1/denp;

%-----Coeficientes R y T (intensidad)-------

Rp=(abs(rp))^2;

Tp=(ql/q1)*(abs(tp))^2;

%------ Matrices ------

GMs=Ms;

GMp=Mp;

%----- Code 2, fitting function------------

function [parameters]=GNTMMs(data_points,nvalues,evalues,lambda,parameters,h,max_counts,min_step)

counts=0;

found=false;

angle=data_points(:,1);

Rsvalues=data_points(:,2);

while (found==false & counts<max_counts)

%------Calculo de los coeficientes para estos parametros------

for k=1:length(angle)

phi0=angle(k);

[~,~,~,Rs,~,Rp,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);

Rfp(k)=Rp;

Rfs(k)=Rs;

%Elementos del jacobiano---

[~,~,~,Rs1p,~,Rp1p,~]=bwTMM([1,parameters(1)+h(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);

[~,~,~,Rs1m,~,Rp1m,~]=bwTMM([1,parameters(1)-h(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);

Rfp1p(k)=Rp1p;

Rfs1p(k)=Rs1p;

Rfp1m(k)=Rp1m;

Rfs1m(k)=Rs1m;

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

[~,~,~,Rs2p,~,Rp2p,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2)+h(2),parameters(3),evalues(1)],phi0,lambda);

[~,~,~,Rs2m,~,Rp2m,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2)-h(2),parameters(3),evalues(1)],phi0,lambda);

Rfp2p(k)=Rp2p;

Rfs2p(k)=Rs2p;

Rfp2m(k)=Rp2m;

Rfs2m(k)=Rs2m;

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

[~,~,~,Rs3p,~,Rp3p,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3)+h(2),evalues(1)],phi0,lambda);

[~,~,~,Rs3m,~,Rp3m,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3)-h(2),evalues(1)],phi0,lambda);

Rfp3p(k)=Rp3p;

Rfs3p(k)=Rs3p;

Rfp3m(k)=Rp3m;

Rfs3m(k)=Rs3m;

end

%------Calculo jacobiano

for r=1:length(Rsvalues)

J(r,1)=(Rfs1p(r)-Rfs1m(r))/(2*h(1));

J(r,2)=(Rfs2p(r)-Rfs2m(r))/(2*h(2));

J(r,3)=(Rfs3p(r)-Rfs3m(r))/(2*h(2));

end

residuals=(Rsvalues-transpose(Rfs))*1000000;

Jt=transpose(J);

stepgn=(1000000*(Jt*J))\(-Jt*residuals);

parameters=parameters+transpose(stepgn)/1000000;

if norm(stepgn)<min_step

found=true;

end

counts=counts+1;

end

end

%------ Code 3--------

close all

clear

clc

data=load('4mM0cs1data.txt');

angle=data(:,1);

Exrp=data(:,2);

Exrs=data(:,3);

Exratio=data(:,4);

%-----known values-----

nvidrio=1.5116+1i*(9.1316)*10^(-9);

evidrio=(1.25)*10^(-6);%nm

nFTO=1.8;

%----------

[parameters]=GNTMMs([angle,Exrs],[nFTO,nvidrio],[evidrio],760,[3,500,200],[0.4,200],50,0.01);

[Torsten]

Why don't you use one of the MATLAB routines for the least squares part, e.g. lsqnonlin ?