Question: extracting data from a contour plot

restart;

with(plots):

# Aim is to understand the general shape of the stability problem, especially for very low and high wavenumbers.

# Parameter values

parvals:= {variphi=0.38, R1=0.0009, R2=8.75, E1=1, E2=1, R4=177.6, C=1.96, k=4, H=7, Pec_i=-120, Phi=1, Q=1.9, R=20};

# Identifying the possible steady positions of the front

frontlocationeq:= (Q*(R*S+1-S)*(H+1/(1-exp(Pec_i)))*exp(C*Pec_i*S/(k*(R*S+1-S)))-H)*(1-exp(Pec_i*(1-S)/(R*S+1-S)))=1;

Svalues:= Student:-Calculus1:-Roots(subs(parvals,frontlocationeq),S);

S01:= op(1,Svalues); S02:= op(2,Svalues); S03:= op(3,Svalues);

# Setting up the dispersion equation

F:=1/(R*S0+1-S0):

w_0:=1/k*Pec_i*F*C:

w_1:=Pec_i*F:

M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i)))):

M_2:=M_1*exp(w_0*S0):

M_3:=w_1/(1-exp(w_1*(S0-1))):

M_4:=M_3*exp(w_1*(S0-1)):

# The dispersion equation is now A = 0, where

A:= -variphi*sigma-2*(((1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(-1/2*(R1-1)*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*M_1*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*C*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*(R1-1)*l*R3*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*exp(S0*(w_0+l))+(-1/2*H*k*l*w_0*R3*(R1-1)*cosh(S0*l)+C*(k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*sinh(S0*l))*(-l*w_0+E1*sigma))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*(R1-1)*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*sigma*variphi*exp(S0*(w_0+l))+cosh(S0*l)*((k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*C+1/2*H*k*sigma*w_0*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*C*k*l*R3*M_1*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+(-l*w_0+E1*sigma)*(-H*k*l*cosh(S0*l)*R3*(R1-1)+C*(M_2*w_0*k+1/2*M_4*w_1)*sinh(S0*l)))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*C*k*sigma*M_1*variphi*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+cosh(S0*l)*((M_2*w_0*k+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(sigma*E2*R1*k+l*C*w_1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(l*w_0+E1*sigma)*(-l*w_0+E1*sigma)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((-C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*exp(w_1*(S0-1)+S0*l)+2*M_4*(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*R1*R2*M_3*sigma*exp(1/2*w_1*(S0-1)+l)*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)))/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/(-w_0^2*l^2+E1^2*sigma^2)/sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0):

# To make progress with marginal stability, assume that the exchange of stabilities holds. The condition for marginal stability is now A0 = 0, where

A0:= subs(sigma=0,A):

# Plot up some cases and see what they look like.

Digits:= 40: # High value of Digits required to avoid numerical noise

# A0 as a function of R3 and l, for the lowest front position

display([contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=black,contours=[0]), contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=blue)]);