How to get tau for every c and plot a c vs tau graph

How to get tau for every c from the equation "tau=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega)" where A, B, C, D1, omega, x, y depends on c.

 clc
  clear 
  close all
  r=3.3;
  K=898;
  alpha=0.045;
  d=1.06;
  h=0.0437;
  theta=0.215;
  for c=.1:.1:.9
       x=d/(alpha*(theta - d*h)*(1 - c));
       y=r*(K-x)*(1+alpha*h*x*(1-c))/(1 - c)*K*alpha;
       A=r*d*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/(alpha*theta*K*(1-c))-d;
       B=d;
       C=-r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/alpha*theta*K*(1-c);
       D1=r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/x*theta*K*(1-c)+r*d^2*(K-x)/alpha*K*theta*(1-c)*x;
       omega=sqrt(((A^2-B^2-2*C)+sqrt((A^2-B^2-2*C)^2-4*(C^2-D1^2)))/2);
       tau=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega);
       plot(c,tau)
     end

Answer 1

[Iman Ansari]

Hi

clc
  clear 
  close all
  r=3.3;
  K=898;
  alpha=0.045;
  d=1.06;
  h=0.0437;
  theta=0.215;
  i=1;
  for c=.1:.1:.9
       x=d/(alpha*(theta - d*h)*(1 - c));
       y=r*(K-x)*(1+alpha*h*x*(1-c))/(1 - c)*K*alpha;
       A=r*d*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/(alpha*theta*K*(1-c))-d;
       B=d;
       C=-r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/alpha*theta*K*(1-c);
       D1=r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/x*theta*K*(1-c)+r*d^2*(K-x)/alpha*K*theta*(1-c)*x;
       omega=sqrt(((A^2-B^2-2*C)+sqrt((A^2-B^2-2*C)^2-4*(C^2-D1^2)))/2);
       tau(i)=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega);
       i=i+1;
  end
plot(.1:.1:.9,tau)

[Atom]

Thank you Ansari. Would you please suggest me why the above programme does not produce any plot.

[Iman Ansari]

In each loop it plot a dot, but you don't use hold on so the graph reset each loop. If using hold on and marker to better see the result:

clc
  clear 
  close all
  r=3.3;
  K=898;
  alpha=0.045;
  d=1.06;
  h=0.0437;
  theta=0.215;
  for c=.1:.1:.9
       x=d/(alpha*(theta - d*h)*(1 - c));
       y=r*(K-x)*(1+alpha*h*x*(1-c))/(1 - c)*K*alpha;
       A=r*d*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/(alpha*theta*K*(1-c))-d;
       B=d;
       C=-r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/alpha*theta*K*(1-c);
       D1=r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/x*theta*K*(1-c)+r*d^2*(K-x)/alpha*K*theta*(1-c)*x;
       omega=sqrt(((A^2-B^2-2*C)+sqrt((A^2-B^2-2*C)^2-4*(C^2-D1^2)))/2);
       tau=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega);
       plot(c,tau,'MarkerSize',15,'Marker','.')
       hold on
     end

[Atom]

Thank you very much. I have ignored that. Thank you again. Mr Ansari, would you please suggest me how to get a continuous line instead of getting dotted plots. Note that I have modified the code as per your suggestion as

 clc
  clear 
  close all
  r=3.3;
  K=898;
  alpha=0.045;
  d=1.06;
  h=0.0437;
  theta=0.215;
  for c=.1:.001:.9
       x=d/(alpha*(theta - d*h)*(1 - c));
       y=r*(K-x)*(1+alpha*h*x*(1-c))/(1 - c)*K*alpha;
       A=r*d*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/(alpha*theta*K*(1-c))-d;
       B=d;
       C=-r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/alpha*theta*K*(1-c);
       D1=r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/x*theta*K*(1-c)+r*d^2*(K-x)/alpha*K*theta*(1-c)*x;
       omega=sqrt(((A^2-B^2-2*C)+sqrt((A^2-B^2-2*C)^2-4*(C^2-D1^2)))/2);
       tau=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega);
       plot(c,tau,'MarkerSize',15,'Marker','.')
       hold on
     end
xlabel('c');
ylabel('tau');

[Iman Ansari]

To plot continuous you need to have two value in each loop for plotting a line. I used tau_0 and c_0 to store previous position in each loop:

clc
  clear 
  close all
  r=3.3;
  K=898;
  alpha=0.045;
  d=1.06;
  h=0.0437;
  theta=0.215;
  for c=.1:.001:.9
       x=d/(alpha*(theta - d*h)*(1 - c));
       y=r*(K-x)*(1+alpha*h*x*(1-c))/(1 - c)*K*alpha;
       A=r*d*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/(alpha*theta*K*(1-c))-d;
       B=d;
       C=-r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/alpha*theta*K*(1-c);
       D1=r*d^2*(alpha*K*h*(1-c)-1-2*alpha*h*(1-c)*x)/x*theta*K*(1-c)+r*d^2*(K-x)/alpha*K*theta*(1-c)*x;
       omega=sqrt(((A^2-B^2-2*C)+sqrt((A^2-B^2-2*C)^2-4*(C^2-D1^2)))/2);
       tau=acos(((omega^2*(D1-A*B)-C*D1)/(B^2*omega^2+D1^2))/omega);
       if c>0.1
           plot([c_0 c],[tau_0 tau])
       end
       tau_0=tau;
       c_0=c;
       hold on
     end
xlabel('c');
ylabel('tau');

[Atom]

Thank you Mr Ansari. Your work helps me a lot. Thank you. Hope I will get help from you in future.