Can anyone tell how can i modify my function as it takes ages to work ?

This code is used to calculate some parameters from an interference pattern image. but acually it takes so long to calculate, maybe more than one hour. here is the code, if there anyone he could help me please.

hamo  = inputdlg({'Enter the function describing the outer shape of the fiber f(x,y)=0','Enter the function describing the shift','Enter the value of the interfringe spacing in micrometre','Enter the value of the used light wavelength in micrometre','Enter the value of the liquid refractive index'},' ',1,{'x^2+y^2-55.01^2','-(0.26743328 -0.13474712*cos(0.057107016*x) +3.7516074e-005*sin(0.057107016*x) -0.066680101*cos(0.11421403*x) -3.8116065e-005*sin(0.11421403*x) -0.033512692*cos(0.17132105*x) -1.2997793e-005*sin(0.17132105*x) -0.018337484*cos(0.22842806*x) +9.8609565e-005*sin(0.22842806*x) -0.0086647943*cos(0.28553508*x) -0.00017735617*sin(0.28553508*x) -0.0034414322*cos(0.3426421*x) +0.00020521148*sin(0.3426421*x) -0.00019872566*cos(0.39974911*x) -0.00015907215*sin(0.39974911*x) +0.0015646411*cos(0.45685613*x) +4.9858462e-005*sin(0.45685613*x))','1','0.6328','1.4996'});
hamam  = inputdlg({'Enter the number of layers you want to divide the fiber into','Enter the initial value of refractive index ni','Enter the accuracy of the refractive index profile (e.g. 0.001)','Enter the value of maximum x-axis (Xm) in micrometre','put 1 for 2 beam or 2 for multiple beam(t)'},' ',1,{'25','1.499','0.0002','110.02','1'});

nn = eval(hamam{1});
ni = eval(hamam{2});
dn = eval(hamam{3});
xm = eval(hamam{4});
dx = xm/(2*nn);
n0 = eval(hamo{5});
h = eval(hamo{3});
lamda = eval(hamo{4});
t = eval(hamam{5});
syms x y;
zz = hamo{2};
p= hamo{1};
f{1} = eval(p);
fx{1} = diff(f{1},'x');
fy{1} = diff(f{1},'y');
for i=2:nn
x= x+x/(nn-i+1);
y= y+y/(nn-i+1);
f{i} = eval(p);
fx{i}=diff(f{i},'x');
fy{i}=diff(f{i},'y');
%pretty(f{i})
%pretty(fx{i})
%pretty(fy{i})
end

ww = xm/2-dx/2;
x =xm -0.5*dx;
n(1) = ni;
for i=1:nn

      RHS = 1010;
      z = eval(zz);

      LHS=(z*lamda)/(1*h);

      while (RHS - LHS)>0.001

      n(i) = n(i)-dn;
      gama = 0;

      [x,y]= solve('x=w',(f{1}));
      w=ww;
      [m,g] = size(x);

      for j = 1:m
          aa(j)= eval(x(j));
          bb(j)= eval(y(j));      
      end
      %The following are the values of (a1,b1)
      a(1) = aa(1);
      b(1) = min(bb);
      bbb(1)=max(bb);
      x = a(1);
      y = b(1);
      c = eval(-1*fy{1}/fx{1});  %here we calculate the constant c1 by the value of f1x
      cetaa =acos(c/(c^2+1)^0.5);
      ceta = (cetaa/pi)*180;
      %dd=cetaa*180/pi

      r = asin((n0/n(1))*sin(cetaa));
      %rr=r*180/pi
      B(1) = cetaa-r;    %Here we calculate Beta = ceta - r.
      gama = gama +B(1);
      if gama > 0
          s(1) = -1*tan(pi/2-gama);
      elseif gama < 0
      s(1) = tan(pi/2+gama);
      end
      t(1) = b(1) - a(1)*s(1);

     if i==1

         if gama ==0
             [x,y]= solve('x-a(1)',(f{1}));
         else
         [x,y]= solve('y-s(1)*x-t(1)',(f{1}));
         end

         [m,g] = size(x);

         for j = 1:m
             aa(j)= eval(x(j));
             bb(j)= eval(y(j));     
         end

         [b(2),I] = max(bb);

         a(2) = aa(I);

      else

              for L = 1:i-1

                  if gama == 0
                      [x,y]= solve('x=a(L)',(f{L+1}));
                  else
                      [x,y]= solve('y-s(L)*x-t(L)',(f{L+1}));
                  end

                  [m,g] = size(x);
                  for j = 1:m
                      aa(j)= eval(x(j));
                      bb(j)= eval(y(j));      
                  end
                  [b(L+1),I] = min(bb);

                  a(L+1) = aa(I);

                  x=a(L+1);
                  y=b(L+1);
                  c = eval(-1*fy{L+1}/fx{L+1});  %here we calculate the constant c1 by the value of f1x
                 cetaa = acos((c*s(L)+1)/(((c^2+1)^0.5)*((s(L)^2+1)^0.5)));
                 if cetaa>pi/2
                      cetaa = pi-cetaa;
                  end
                 %ee=cetaa*180/pi

                  r = asin((n(L)/n(L+1))*sin(cetaa));
                  B(L+1) = cetaa-r;    %Here we calculate Beta = ceta - r.
                  gama = gama +B(L+1);
                  if gama > 0
                  s(L+1) = -1*tan(pi/2-gama);
                  elseif gama<0
                  s(L+1) = tan(pi/2+gama);
                  end
                  t(L+1) = b(L+1) - a(L+1)*s(L+1);

              end

              for L = i:(2*i-1)

                  if gama == 0
                      [x,y]= solve('x=a(L)',(f{2*i-L}));
                  else
                      [x,y]= solve('y-s(L)*x-t(L)',(f{2*i-L}));
                  end

                  [m,g] = size(x);
                  for j = 1:m
                      aa(j)= eval(x(j));
                      bb(j)= eval(y(j));      
                  end

                  [b(L+1),I] = max(bb);

                  a(L+1) = aa(I);

                  if L<(2*i-1)
                  x=a(L+1);
                  y = b(L+1);
                  c = eval(-1*fy{2*i-L}/fx{2*i-L});  %here we calculate the constant c1 by the value of f1x
                  cetaa = acos((c*s(L)+1)/(((c^2+1)^0.5)*((s(L)^2+1)^0.5)));
                  if cetaa>pi/2
                      cetaa = pi-cetaa;
                  end
                  %gg = cetaa*180/pi
                  r = asin((n(2*i-L)/n(2*i-L-1))*sin(cetaa));
                  B(L+1) = cetaa-r;    %Here we calculate Beta = ceta - r.
                  gama = gama +B(L+1);
                  if gama>0
                  s(L+1) = -1*tan(pi/2-gama);
                  elseif gama<0
                   s(L+1) = tan(pi/2+gama);
                  end
                  t(L+1) = b(L+1) - a(L+1)*s(L+1);
                  end

              end

      end

       L0 = abs(b(2*i)-b(1));
       LLi =((a(i+1)-a(i))^2+(b(i+1)-b(i))^2)^0.5;
       ssss(i)=LLi;
       L=0;
  %      llllll=0;
       for q=1:i-1
           LL(q)= ((a(q+1)-a(q))^2+(b(q+1)-b(q))^2)^0.5 + ((a(2*i-q+1)-a(2*i-q))^2+(b(2*i-q+1)-b(2*i-q))^2)^0.5;
           L = L + LL(q)*n(q);
           llllll(q)=((a(q+1)-a(q))^2+(b(q+1)-b(q))^2)^0.5;             %here we calculate lllll  which is the real length
           lllllll(q)=((a(2*i-q+1)-a(2*i-q))^2+(b(2*i-q+1)-b(2*i-q))^2)^0.5;
  end

       RHS = n(i)*LLi+L-n0*L0;

      end
      n(i+1)=n(i)+0.004
  % n(i+1)=n(i)
      pol(i)=(3*(n(i)^2-1))/(4*pi*(n(i)^2+2));
      nm(i)=n(i);
      x=xm/2+a(2*i);
      ww = ww-dx;

end    

%  navg1=n(1);
%      for i=2:nn
%     navg1= navg1+n(i);
%      end
% navg=navg1/nn;
% 
% p1=(3*(navg^2-1))/(4*pi*(navg^2+2));
[file,path] = uiputfile('*.txt','Save As');
df = fullfile(path,file);
mm = fopen(df, 'wt');

fprintf(mm, '%s                 %s\n','x','n');
for i=1:nn
x = -(nn-0.5*(2*i-1))*dx;
a = n(i);
fprintf(mm, '%f         %f\n',x,a);
end
for i=1:nn
x = (nn-0.5*(2*i-1))*dx;
a = n(i);
fprintf(mm, '%f         %f\n',x,a);
end

fclose(mm)

[file,path] = uiputfile('*.txt','Save As');
df = fullfile(path,file);
mm = fopen(df, 'wt');

fprintf(mm, '%s                 %s\n','x','l');
for i=1:nn-1
x = -i;
a = llllll(i);
fprintf(mm, '%f         %f\n',x,a);
end
% if i==nn
x = nn;
a = LLi/2;
fprintf(mm, '%f         %f\n',x,a);
x = nn;
a = LLi/2;
fprintf(mm, '%f         %f\n',x,a);
% end
for i=1:nn-1
x = i;
a = lllllll(i);
fprintf(mm, '%f         %f\n',x,a);
end

fclose(mm)

[Stephen23]

Use the symbolic tools (e.g. subs, simplify, double, etc) rather than eval: it is the wrong tool for the job!

https://www.mathworks.com/matlabcentral/answers/228557-experts-of-matlab-how-did-you-learn-any-advice-for-beginner-intermediate-users

Answer 1

[Walter Roberson]

You should avoid using eval() as it prevents the code from being optimized. You are also using eval() with symbolic expressions (it looks like) which is a problem because the language of symbolic expressions is not exactly MATLAB.

As you are using R2014b (thanks for filling in that information!) you should use sym() to convert character vectors into symbolic expressions, and you should use str2fun() to convert character vectors into anonymous functions.

syms x y;

You do that after you have first eval()'d the hamo entries, so the x and y values in those entries is going to be evaluated using whatever value of x and y happens to be in the workspace. Then you assign symbolic values to x and y, confusing people who might have thought they knew what x and y meant. And you have a bunch of lines like

[x,y]= solve('x-a(1)',(f{1}))

which assigns new values again to x and y. Note that when you do that, then any symbolic expression that was created with symbolic x or y previously does not automatically get updated, so you can end up with symbolic x in expressions that do not refer to the x that resulted from the solve(). This is confusing and often results in errors.

My advice is to never reassign the value of a symbolic variable for a completely different use once you have used the value of the variable -- not unless you are going to go ahead and overwrite all of the places you used the value before.