function retArray = mobCal5Si(fileToOpen,VgIndex,VdIndex,IdIndex,k_max_in,l_w_in)
%arguments- filename,gate voltage index,drain voltage index,drain current index,double sweep,length_width 
%written for calulating the mobility
%first argument is the name of the file ( usally a CSV got from DC probe
%station
%header is the number of header lines, before the actual data starts holds
%information about the number of data rows the file has
%VgIndex is in the data matrix which one is Vg ( 1-4) other Indexes are
%similar to VG
%Fits a singlular value of mobility and N_0

header = returnHeaderLines(fileToOpen,'DataName');
inData=importdata(fileToOpen,',',header);
data = inData.data;
intext = inData.textdata;

%%% search for hysteresis
if ~isempty(cell2mat(regexp(intext(16),'Double', 'once')))
    k_max = 2;
elseif ~isempty(cell2mat(regexp(intext(16),'Single', 'once')))
    k_max=1;
else
    disp('wrongly read/ mialigned and using passed Value')
    k_max = k_max_in;
end

%% getting the dimensions of the sweep from the CSV file
s2 = regexp(intext(header-1,1),',','split');
d2 = str2double(s2{1}{2});
s1 = regexp(intext(header-2,1),',','split');
d1_in = str2double(s1{1}{2});
if (length(inData.data) ~= d1_in*d2)
    disp('truncated file');
    return;
end
mycolor =distinguishable_colors(d2+2); 
%xaxis = 1:1:d1;
% plot(data(:,1),data(:,3));
d1 = floor(d1_in/k_max);
%for z=1:1:2
%       data(:,:,1) = data_in(1:d1,:);
%       data(:,:,2) = data_in(d1+1:d1_in,:);
Vg = zeros(d1,d2,k_max);
Id = zeros(d1,d2,k_max);
Vd = zeros(d1,d2,k_max);
R = zeros(d1,d2,k_max);
R_fit = zeros(d1,d2,k_max);
g = zeros(d1,d2,k_max);
n_org = zeros(d1,d2,k_max);
n = zeros(d1,d2,k_max);
mu = zeros(d1,d2,k_max);
mu3 = zeros(d1,d2,k_max);
mu2 = zeros(d1,d2,k_max);
maxres = zeros(d2,k_max);
index = zeros(d2,k_max);
maxResVg = zeros(d2,k_max);


for j=1:1:d2
    for k = 1:1:k_max
        for i=1:1:d1
            Vg(i,j,k)=data((j+k-2)*d1 + i,VgIndex);
            Vd(i,j,k)= data((j+k-2)*d1+ i,VdIndex);
            Id(i,j,k) = data((j+k-2)*d1 + i,IdIndex);
            R(i,j,k) = abs(Vd(i,j,k)/Id(i,j,k));
            if  i == 1 
                g(i,j,k) = Id(i,j)/Vg(i,j);
            else
                g(i,j,k) = (Id(i,j,k) - Id(i-1,j,k))/(Vg(i,j,k) - Vg(i-1,j,k));
            end
        end
[maxres(j,k),index(j,k)] = max(R(:,j,k));
maxResVg(j,k) = data((j+k-2)*d1 + index(j,k),VgIndex);
    end
end

%for i = 1:1:d2
%    plot(data((d1*(i-1) + 1):d1*i,1),abs(R((d1*(i-1)+1):d1*i)),cstring(i));
%end

% calculation of surface carrier density
ep_air = 1;
ep_ox = 3.9;
%ep_ox = 9;% for AlGaN
t_air = 0;
ec=1.6e-19;
ep_0 = 8.854E-14;
t_ox = 300e-7;% 30 for AlGAN, and 300 for SiO2
c_air = ep_0 * ep_air/t_air; %%%% F/m2
c_ox = ep_0 * ep_ox/t_ox; %% F/m2
l_w = l_w_in;
h = 6.6260E-34;
cond_min = 4 * (ec^2)/(pi*h);
if (t_air ~= 0 && t_ox ~=0 ) 
c_eff = c_air*c_ox/(c_air + c_ox); %%F/m2
elseif (t_air == 0)
    c_eff = c_ox;
else
    c_eff = c_air;
end

if(k_max>1)
sse_hist=0;
R_hist=zeros(d2);
for j=1:1:d2
    for i = 1:1:d1
        R_hist(i,j) = (R(i,j,1) - R(d1-i+1,j,2))^2;
        sse_hist = sse_hist + R_hist(i,j);
    end
   hist_trap_density(j) = ((maxResVg(j,2) - maxResVg(j,1))*c_eff)/2*ec; 
end
rmse_hist_inter(j) = sqrt(mean(R_hist(:,j)));
rmse_hist = mean(rmse_hist_inter);
mean_trap = mean(hist_trap_density);
else
    sse_hist = 0;
    mean_trap = 0;
    rmse_hist = 0;
end

Rc = zeros(d2,2);

mu_fit = zeros(d2,2);

n0 = zeros(d2,2);

gof_r = zeros(d2,2);

gof_rmse = zeros(d2,2);
gof_sse = zeros(d2,2);


maxMob = zeros(d2,2);

for j=1:1:d2
    for k=1:1:k_max
        for i=1:1:d1
        n_org(i,j,k) = (Vg(i,j,k)*c_eff)/ec;
        n(i,j,k) =(Vg(i,j,k) - maxResVg(j,k))*c_eff/ec;
        mu3(i,j,k) = (l_w)* g(i,j,k)/(c_eff*Vd(i,j,k));   %% max transconductance model
        mu2(i,j,k) = (l_w* g(i,j,k))/(n(i,j,k)*ec);
        mu(i,j,k) = l_w/(R(i,j,k) * n(i,j,k) * ec); %% drude model
        end;
    r_vec_n = R(:,j,k);
    n_vec_n = n(:,j,k);

% Fit this model using new data
%ntype = fit(n_vec_n,r_vec_n,ft_,fo_);
    try
    [ntype,gof_n] = fitMobFun(n_vec_n,r_vec_n);
    R_fit(:,j,k) = feval(ntype,n_vec_n);
    Rc(j,k)= ntype.c/2;
    mu_fit(j,k) = (ntype.a)/l_w;
    n0(j,k) = ntype.b;
    gof_r(j,k) = gof_n.rsquare;
    gof_rmse(j,k) = gof_n.rmse;
    gof_sse(j,k) = gof_n.sse;
    catch err
    disp(strcat('error in fitting ntype',fileToOpen));
    disp(err);
    Rc(j,k) = 0;
    mu_fit(j,k) = 0;
    n0(j,k)=0;
    gof_r(j,k) = 0;
    gof_rmse(j,k) = 0;
    gof_sse(j,k) = 0;
    end
    if ~(gof_n.rsquare > 0.5 && gof_n.rsquare <1) 
        Rc(j,k) = 0;
        mu_fit(j,k) = 0;
        n0(j,k)=0;
    end
%%%% to plot the fitted curved
%ntype_fit2 = feval(ntype,n_vec_n);
%

    

    maxMob(j,k) = max(mu(mu(:,j,k)<inf,j,k));

    end
end

% h=figure;
% hold on;
% legendString = cell(1,d2);
% for i = 2:1:d2
%     for k=1:1:k_max
%     plot(Vg(:,i,k),abs(mu3(:,i,k)),'Color',mycolor(i,:),'LineWidth',3);
%     legendString{i} = sprintf('Vd = %f',Vd(1,i));
%     end
% end
% xlabel('gate voltage');
% ylabel('Mobility');
% title(fileToOpen);
% legend(legendString{1:d2});

% hold off;
% print(h,'-depsc',strcat(fileToOpen,'-maxT.eps'));
% close(h);
h=figure;
hold on;
%legendString = cell(1,d2);
for i = 2:1:d2
    for k =1:1:k_max
        subplot(2,2,1)
        hold on
        plot(Vg(:,i,k),abs(mu3(:,i,k)),'Color',mycolor(i,:),'LineWidth',2);
        xlabel('gate voltage');
        ylabel('Mobility');
        subplot(2,2,2)
        hold on
        plot(Vg(:,i,k),abs(mu(:,i,k)),'Color',mycolor(i,:),'LineWidth',2);
        xlabel('gate voltage');
        ylabel('Mobility');
        subplot(2,2,3)
        hold on
        plot(Vg(:,i,k),R(:,i,k),'Color',mycolor(i,:),'LineWidth',2);
        xlabel('gateVoltage');
        ylabel('Resistance');
        subplot(2,2,4)
        hold on
        plot(Vg(:,i,k),R_fit(:,i,k),'Color',mycolor(i,:),'LineWidth',2);
        xlabel('gate voltage');
        ylabel('Fitted Resistance');
             
    end
end

%legend(legendString{1:d2},'Location','NorthEastOutside');
hold off;
%print(h,'-depsc',strcat(fileToOpen,'.eps'));
print(h,'-djpeg',strcat(fileToOpen,'.jpeg'));
close(h);

minCond = zeros(1);
Dp = zeros(1);
maxMobMean = zeros(1);
meanRc = zeros(1);

meanMu = zeros(1);

meanN0 = zeros(1);

meanGof_r = zeros(1);

meanGof_rmse = zeros(1);

meanGof_sse = zeros(1);


fp2 = fopen(strcat(fileToOpen,'_allInfo.dat'),'a');
fprintf(fp2,'\nName of the file-%s\t executed on %s\n',fileToOpen,datestr(now));
for k =1:1:k_max
minCond(k) = 1/(mean(maxres(:,k))*cond_min);
Dp(k) = mean(maxResVg(:,k));
maxMobMean(k) = mean(maxMob(:,k));
meanRc(k) = mean(Rc(:,k));

meanMu(k) = mean(mu_fit(:,k));
meanN0(k) = mean(n0(:,k));

meanGof_r(k) = mean(gof_r(:,k));

meanGof_rmse(k) = mean(gof_rmse(:,k));

meanGof_sse(k) = mean(gof_sse(:,k));


fprintf(fp2,'\nmean Rc--%f\tminCond--%f\tDiracPointPosition--%f\nmeanMu--%f\tmaxMobility--%f\nmeanN0--%e\nmeanGof_r--%f\tmeanGof_rmse--%f\tmeanGof_sse--%e\nsse_hist--%e\trmse_hist--%e\t mean_trap_density --%e\n\n',meanRc(k),minCond(k),Dp(k),meanMu(k),maxMobMean(k),meanN0(k),meanGof_r(k),meanGof_rmse(k),meanGof_sse(k),sse_hist,rmse_hist,mean_trap);

end




fclose(fp2);
retArray = [meanRc,meanMu,meanN0,maxMobMean,minCond,Dp,meanGof_r,meanGof_rmse,meanGof_sse,sse_hist,rmse_hist];

return;
