%Parameterized seaching method, including iterating beta%% Picks an initial phi, runs collocation to get a profile and% uses that to calc eta.% Then finds a point (phi,eta) along the curve that is some% predetermined DISTANCE from the previous point% The variable psi is the angle of the step along a curve with respect% to the x-axis% The radius of search rE is on the log scale, it is the radius of search% on the actual figure, which is then convered to linear scale.% Written by Brad Schwab, UW, 5/2018.% Slight editing by jbr, 5/2018.%% Eliminate need for global variables% Struct p contains the following pararmeters:% R A B phi_prev phi_ini eta_prev gamma beta cS rE I%tic%Changeable parameters for optimizationp.gamma=30;n=60;%number of collocation pointsnpoint=75;%max number of steps along the curvebetavals=[0.6;0.4;0.3;0.2;0.1;0;-0.8];rE0=.15*ones(length(betavals),1);rE0(5)=.1;%Search will skip the top of a curve at rE0=.15 for beta=0%Shape factor and boundary conditionRad=3;p.cS=1;%Adds aysmptotic line and verticle line seen in book Figure 7.19xline1=[.00110];yline1=[1000.1];xline2=[.01.01];yline2=[.2500];%Suppress output from fsolveoptions=optimset('Display','off');%Get collocation points, note the same n is used for all beta/phi values[R,A,B,Q]=colloc(n-2,'left','right');p.R=R*Rad;p.A=A/Rad;p.B=B/(Rad*Rad);p.Q=Q*Rad;%Storage of (phi,eta) values and initializes first runnbeta=numel(betavals);phistore=cell(nbeta,1);etastore=cell(nbeta,1);p.phi_ini=5E-4;para_end=.01;%Initial guesses for fsolveIn01=ones(n,1);%need an inital eta to start parametric searchloglog(xline1,yline1);holdonloglog(xline2,yline2);forj=1:nbetap.beta=betavals(j);p.I=calcI(p.gamma,p.beta);p.rE=rE0(j);%initial radius is set to base, changed in loop%Phi is fixed at phi_ini for initial searchp.phi_prev=p.phi_ini;%start searching at phi = .01, changed in looppsi_prev=0;delpsi=0;In0=[ones(n,1);0];fori=1:npointif(i==1)%First need to find an eta for the initial phi[ini,resid,info]=fsolve(@(x)eta_initial(x,p),In01,options);deriv=p.A*ini;last_der=deriv(end);p.eta_ini=last_der/(p.phi_ini*p.I)^2;%Set storage variables for first searchp.eta_prev=p.eta_ini;phistore{j}(i)=p.phi_ini;etastore{j}(i)=p.eta_ini;end%if%Next the parametric search starts[x,resid,info]=fsolve(@(x)eta_calc(x,p),In0,options);psi=x(end);eta=10^(p.rE*sin(psi)+log10(p.eta_prev));phi=10^(p.rE*cos(psi)+log10(p.phi_prev));In0=x;%Next set of c are the final values of last setIn0(end)+=delpsi;%Next psi is calculated as a fucntion%of the last two psi values%Set storage values for this iphistore{j}(i)=phi;etastore{j}(i)=eta;p.phi_prev=phi;p.eta_prev=eta;delpsi=psi-psi_prev;psi_prev=psi;%Dynamically change radius of search based on change in psi of%last seach rE in (0,rE0)p.rE=rE0(j)*(1-abs((delpsi)/pi))^20;if(abs(log10(p.eta_prev*p.phi_prev))<para_end)breakend%ifend%for%Graph curve for this betaloglog(phistore{j},etastore{j},'o')table{j}=[phistore{j};etastore{j}]';end%fortable2=[xline1(:),yline1(:),xline2(:),yline2(:)];%Labels on main figurexlabel('phi')ylabel('eta')axis([1e-4,10,0.1,1e3])%tocsavewh.dattabletable2