i need to understand this code

function result=bend_extend_twist_function(x,psi,theta,lengths,numSegsInit,numBendSegsInit,numSteps,Ri,Rm,Ro,e,area,mu1,mu2,muBar,phi,segTypes)

pressure=x(1:numSteps);
alpha=x(numSteps+1:numSteps+numSegsInit);
twistSegLengthsInit=x(numSteps+1+numSegsInit:end-numBendSegsInit);
bendSegLengthsInit=x(end-numBendSegsInit+1:end);

%%
options = optimset('display','off','MaxFunEvals',10000,'MaxIter',10000,'tolfun',1e-12,'tolcon',1e-12,'tolx',1e-12);

eps1 = (Rm - Ri)/Ri;
eps2 = (Ro - Rm)/Rm;

Lf = abs(1./sin(alpha)); %length of fiber 
num = max(abs(round(12*sin(alpha))),1); %number of fibers

Vf = num*area.*Lf; %fiber volume
Ve_outer = pi*(Ro^2-Rm^2); %elastomer volume in outer layer

t1=Rm-Ri; % actuator inner layer wall thickness
t2=Ro-Ri; % actuator total wall thickness

numSegs=length(alpha); %number of fiber angles ie segments 

%% Calculate initial segment lengths, and required lengths of extending segments at subsequent pressures

extendCount=0; bendCount=0; twistCount=0;
if segTypes(1)==2 
    twistCount=twistCount+1;
    bendCount=bendCount+1;
    segLengthsInit(1) = twistSegLengthsInit(twistCount);
    segLengthsInit(2) = bendSegLengthsInit(bendCount);

elseif segTypes(1)==0
    extendCount=extendCount+1;
    bendCount=bendCount+1;
    segLengthsInit(1) = lengths(extendCount,1)-bendSegLengthsInit(bendCount)/2;
    segLengths(extendCount,:) = lengths(extendCount,:)-bendSegLengthsInit(bendCount)/2-0.5*Ri*(psi(bendCount,:)*pi/180)*(cos(phi)+1)-0.5*t2*psi(bendCount,:)*pi/180;
    segLengthsInit(2) = bendSegLengthsInit(bendCount);
end

for i=3:length(segTypes)
    if segTypes(i)==0
        extendCount=extendCount+1;
        segLengthsInit(i) = lengths(extendCount,1)-bendSegLengthsInit(bendCount)/2;
        segLengths(extendCount,:) = lengths(extendCount,:)-bendSegLengthsInit(bendCount)/2-0.5*Ri*(psi(bendCount,:)*pi/180)*(cos(phi)+1)-0.5*t2*psi(bendCount,:)*pi/180;

    elseif segTypes(i)==1
        bendCount=bendCount+1;
        segLengthsInit(i) = bendSegLengthsInit(bendCount);
        segLengthsInit(i-1)=segLengthsInit(i-1)-bendSegLengthsInit(bendCount)/2;
        
        segLengths(extendCount,:) =  segLengths(extendCount,:)-bendSegLengthsInit(bendCount)/2-0.5*Ri*(psi(bendCount,:)*pi/180)*(cos(phi)+1)-0.5*t2*psi(bendCount,:)*pi/180;
    end
end


%%
psi_per_L = (psi./repmat(bendSegLengthsInit',1,size(psi,2)))*pi/180; %bend per unit length

R=1./psi_per_L+Ri*cos(phi); %radius of curvature of bending segment

%% With the current set of variables: 
% for each extending segment, find the extension achieved.
% for each bending segment, find the pressure at which the given bend angle is achieved.
% for each twisting segment, find the twist angle achieved

output=zeros(numSegs,numSteps);Lz=zeros(numSegs,numSteps);a=zeros(numSegs,numSteps);tau=zeros(numSegs,numSteps);
counter=0;

for i=1:numSegs %for each segment
  
    if segTypes(i)==0 %if it is an extending segment
        
        Vtot_outer = Ve_outer+Vf(i)+Vf(i); %total volume of outer layer
        rho = (num(i)/Vtot_outer)*e*area*Lf(i)/2;
        rho_elast_iso = mu1;
        rho_elast_aniso = mu1*Ve_outer./Vtot_outer;
        
        %solve at initial pressure:
        P=pressure(1);
        x0=[1,Ri,0.001];
        
        f = @(x)numerics_Taylor_2fam(Ri,Rm,x,P,rho_elast_iso,rho_elast_aniso,alpha(i),rho,rho,eps1,eps2);
        
        [x,fval] = fsolve(f,x0,options);
        
         output(i,1)=x(1)*segLengthsInit(i);%segment length
         Lz(i,1)=x(1);%segment extension
         a(i,1)=x(2);%segment radius
         tau(i,1)=x(3);%segment twist
          
    %solve at subsequent pressures:
    for j=2:numSteps 
        
        P=pressure(j);
        x0=[Lz(i,j-1),a(i,j-1),tau(i,j-1)];
        
        f = @(x)numerics_Taylor_2fam(Ri,Rm,x,P,rho_elast_iso,rho_elast_aniso,alpha(i),rho,rho,eps1,eps2);
        
        [x,fval] = fsolve(f,x0,options);
        
         output(i,j)=x(1)*segLengthsInit(i);
         Lz(i,j)=x(1);
         a(i,j)=x(2);
         tau(i,j)=x(3);

    end
    tau(i,1)=0;
        
    elseif segTypes(i)==1 %if it's a bending segment
        
        counter=counter+1;
        
        Vtot_outer = Ve_outer+Vf(i)+Vf(i);
        
        rho1=(muBar/2)*Ve_outer/Vtot_outer;
        rho2=num(i)*e*area.*Lf(i)/(2*Vtot_outer);
        rho3=(muBar/2)*Ve_outer/Vtot_outer;
        rho4=num(i)*e*area.*Lf(i)/(2*Vtot_outer);
    
        for j=1:numSteps

            Lc=@(x,y)(R(counter,j)-(Ri+x).*cos(y))/(R(counter,j)-Ri*cos(phi));

            I4=@(x,y)cos(alpha(i))^2+Lc(x,y).^2*sin(alpha(i))^2;

            s1=@(x,y)muBar*(Lc(x,y)-1./(Lc(x,y).^3));
            s2=@(x,y)-2*rho1./(Lc(x,y).^3)+2*Lc(x,y)*rho1+2*(-1+sqrt(I4(x,y))).*Lc(x,y)*rho2*sin(alpha(i))^2./sqrt(I4(x,y))+2*(-1+sqrt(I4(x,y))).*Lc(x,y)*rho2*sin(alpha(i))^2./sqrt(I4(x,y));
            s3=@(x,y)muBar*(Lc(x,y)-1./(Lc(x,y).^3));
            s4=@(x,y)-2*rho3./(Lc(x,y).^3)+2*Lc(x,y)*rho3+2*(-1+sqrt(I4(x,y))).*Lc(x,y)*rho4*sin(alpha(i))^2./sqrt(I4(x,y))+2*(-1+sqrt(I4(x,y))).*Lc(x,y)*rho4*sin(alpha(i))^2./sqrt(I4(x,y));

            f1=@(x,y)s1(x,y).*(Ri+x).*(Ri*cos(phi)-(Ri+x).*cos(y));
            f2=@(x,y)s2(x,y).*(Ri+x).*(Ri*cos(phi)-(Ri+x).*cos(y));
            f3b=@(x,y)s3(x,y).*(Ri+x).*(Ri*cos(phi)-(Ri+x).*cos(y));
            f4b=@(x,y)s4(x,y).*(Ri+x).*(Ri*cos(phi)-(Ri+x).*cos(y));

            M1=2*quad2d(f1,0,t1,pi/2,pi);
            M2=2*quad2d(f2,t1,t2,pi/2,pi);
            M3=2*quad2d(f3b,0,t1,0,phi)+2*quad2d(f3b,0,t1,phi,pi/2);%
            M4=2*quad2d(f4b,t1,t2,0,phi)+2*quad2d(f4b,t1,t2,phi,pi/2);%

            Ma=-2*Ri^3*((sin(phi))^3/3+0.5*cos(phi)*(-phi+pi+cos(phi)*sin(phi)))-(Ri^3)*(-phi*cos(phi)+0.75*sin(phi)+(1/12)*sin(3*phi));%negative moment above and below NA
            
            output(i,j)=-(M1+M2+M3+M4)/Ma;
        end
        
    elseif segTypes(i)==2 %if it's a twisting segment
        
        Vtot_outer = Ve_outer+Vf(i);
        
        rho = (num(i)/Vtot_outer)*e*area*Lf(i)/2;
        
        rho_elast_iso = mu1;
        rho_elast_aniso = mu1*Ve_outer./Vtot_outer;
        
        P=pressure(1);
        x0=[1,Ri,0.001];
        
        f = @(x)numerics_Taylor(Ri,Rm,x,P,rho_elast_iso,rho_elast_aniso,alpha(i),rho,eps1,eps2);
        
        [x,fval] = fsolve(f,x0,options);
        
         Lz(i,1)=x(1);
         a(i,1)=x(2);
         tau(i,1)=x(3);
         output(i,1)=abs(tau(i,1)*segLengthsInit(i)*Lz(i,1));
        
        for j=2:numSteps %for each other pressure value

            P=pressure(j);
            x0=[Lz(i,j-1),a(i,j-1),output(i,j-1)];

            f = @(x)numerics_Taylor(Ri,Rm,x,P,rho_elast_iso,rho_elast_aniso,alpha(i),rho,eps1,eps2);

            [x,fval] = fsolve(f,x0,options);

             Lz(i,j)=x(1);
             a(i,j)=x(2);
             tau(i,j)=x(3);
             output(i,j)=abs(tau(i,j))*segLengthsInit(i)*Lz(i,j); %calculate total twist (tau*current length)

        end
                
    else sprintf('error in "type"')
        
    end
end

%% Compare the achieved values to the desired values
count=1;
for i=1:length(segTypes)
    if segTypes(i)==0
      matchThis(i,:) = segLengths(count,:);
      count=count+1;
    elseif segTypes(i)==1
        matchThis(i,:) = pressure;
    elseif segTypes(i)==2
        matchThis(i,:) = theta*pi/180;
    end
end


result=matchThis-output;

%If a segment is a twisting or bending segment, multiply by a weight:
for i=1:length(segTypes)
    if segTypes(i)==1 
        result(i,:) = 1000*result(i,:); 
    elseif segTypes(i)==2
        result(i,:) = 100*result(i,:);
    end
end