how can I find the intersection of two surface

Hi, I wrote the code below to find the intersection of two surfaces "cdot: and "ctdot", but it has an error : Z must be a matrix, not a scalar or vector.

Could you please tell me how solve it, thanks for any help.

vplc=0.16;

delta=0.1;

Ktau=0.045;

Kc=0.1;

K=0.0075;

Kp=0.15;

gamma=5.5;

Kb=0.4;

alpha0=delta*6.81e-6/(0.002);

alpha1=delta*2.27e-5/(0.002);

Ke=7;

Vs=0.002;

Ks=0.1;

Kf=0.18;

kplc=0.055;

ki=2;

tau_max=0.176;

Vss=0.044;

A=(-(Vss.*c.^2)./(Ks.^2))+((Vs.*K.*gamma.^2.*ct.^2)./(Ks.^2))+alpha0+alpha1.*((Ke.^4)./(Ke.^4+(gamma.*ct).^4));

p=(vplc.*c.^2./(c.^2+kplc.^2))./ki;

h=(-(0.4.*A.*((Kc.^4).*(Kp.^2))./((p.^2.*c.^2.*gamma.*ct.*Kf))));

G1=alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4));

G2=((1-h)./tau_max).*c.^4;

Fc=(4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2);

Fct=((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2);

Fh=(gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4);

cdot=@(ct,c)(((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2)).*(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4)))+((gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4)).*(((1-h)./tau_max).*c.^4);

ctdot=@(ct,c)(-(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4))).*((4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2)));

ct = 0:.1:2.5;

c = 0:.1:2.5;

[Ct,C] = meshgrid(ct,c);

figure

surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))

hold on

contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',2)

Answer 1

[Star Strider]

MATLAB is case-sensitive so ‘ct’ ~= ‘Ct’ (and so for the others) in the ‘cdot’ and ‘ctdot’ calls.

I would like to run this to demonstrate that, however ‘c’ is nowhere to be found.

vplc=0.16;
delta=0.1;
Ktau=0.045;
Kc=0.1; 
K=0.0075;
Kp=0.15;
gamma=5.5;
Kb=0.4;
vss=0.044;
alpha0=delta*6.81e-6/(0.002);
alpha1=delta*2.27e-5/(0.002);
Ke=7;
Vs=0.002;
Ks=0.1;
Kf=0.18;
kplc=0.055;
ki=2;
tau_max=0.176;
Vss=0.044;

ct = 0:.1:2.5;
c = 0:.1:2.5;

A=(-(vss.*c.^2)./(ks.^2))+((Vs.*K.*gamma.^2.*ct.^2)./(ks.^2))+alpha0+alpha1.*((Ke.^4)./(Ke.^4+(gamma.*ct).^4));

Unrecognized function or variable 'ks'.

p=(vplc.*c.^2./(c.^2+kplc.^2))./ki;
h=(-(0.4.*A.*((Kc.^4).*(Kp.^2))./((p.^2.*c.^2.*gamma.*ct.*Kf))));
G1=alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4));
G2=((1-h)./tau_max).*c.^4;
Fc=(4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2);
Fct=((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2);
Fh=(gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4);

cdot=@(ct,c)(((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2)).*(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4)))+((gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4)).*(((1-h)./tau_max).*c.^4);
ctdot=@(ct,c)(-(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4))).*((4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2)));

% ct = 0:.1:2.5;
% c = 0:.1:2.5;
[Ct,C] = meshgrid(ct,c);

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',2)  

.

[M]

Thanks for your answer. c is defined like below:

ct = 0:.1:2.5;

c = 0:.1:2.5;

[Ct,C] = meshgrid(ct,c);

May be you have to put the upper code before all of the functions?

I don't have athis problem when I run the code, the only error is one that I wrote in question.

[Star Strider]

O.K. Those were out of sequence in the code, so I put them above where the error occurred.

Now, ‘ks’ is missing, and I do not see it anywhere in the code. (The reason you did not get that error is that the variable was probably defined earlier and I do not have the earlier code. That could also be the reason I do not have ‘ks’.)

I will keep working on this when all the missing variables appear. (I do not know how may others might be missing, since I am discovering them one at a time.)

[M]

@Star Strider you're right...I again check them some of them was written with capital, so I correct them and edit the question code.

please let me know if it has other error because they are not shown me.

[Star Strider]

vplc=0.16;
delta=0.1;
Ktau=0.045;
Kc=0.1; 
K=0.0075;
Kp=0.15;
gamma=5.5;
Kb=0.4;
alpha0=delta*6.81e-6/(0.002);
alpha1=delta*2.27e-5/(0.002);
Ke=7;
Vs=0.002;
Ks=0.1;
Kf=0.18;
kplc=0.055;
ki=2;
tau_max=0.176;
Vss=0.044;

ct = 0:.1:2.5;
c = 0:.1:2.5;

A=(-(Vss.*c.^2)./(Ks.^2))+((Vs.*K.*gamma.^2.*ct.^2)./(Ks.^2))+alpha0+alpha1.*((Ke.^4)./(Ke.^4+(gamma.*ct).^4));
p=(vplc.*c.^2./(c.^2+kplc.^2))./ki;
h=(-(0.4.*A.*((Kc.^4).*(Kp.^2))./((p.^2.*c.^2.*gamma.*ct.*Kf))));
G1=alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4));
G2=((1-h)./tau_max).*c.^4;
Fc=(4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2);
Fct=((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2);
Fh=(gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4);

cdot=@(ct,c)(((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2)).*(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4)))+((gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4)).*(((1-h)./tau_max).*c.^4);
ctdot=@(ct,c)(-(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4))).*((4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2)));

[Ct,C] = meshgrid(ct,c);

figure
surf(Ct,C,cdot(Ct,C))
hold on
surf(Ct,C,zeros(size(C)), 'FaceColor','r', 'FaceAlpha',0.5)
hold off
xlabel('X')
ylabel('Y')
title('cdot')

figure
surf(Ct,C,ctdot(Ct,C))
hold on
surf(Ct,C,zeros(size(C)), 'FaceColor','r', 'FaceAlpha',0.5)
hold off
xlabel('X')
ylabel('Y')
title('ctdot')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
xlabel('X')
ylabel('Y')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
xlabel('X')
ylabel('Y')
ylim([0 0.2])
zlim([-1 1])
title('Zoomed')

Mtx = cdot(Ct,C)-ctdot(Ct,C);
[minMtx,idx] = min(Mtx(:))

minMtx = 1.8648e-05

idx = 651

[row,col] = ind2sub(size(Mtx), idx)

row = 1

col = 26

idx = find(Mtx == min(Mtx(:)))

idx = 651

[row,col] = ind2sub(size(Mtx), idx)

row = 1

col = 26

I am not certain that the surfaces actually intersect. If so, they do so in a very narrow range of ‘Y’ although likely for all ‘X’ in that region. (I plotted a plane equal to zero in both plots to see if there could be an intersection. Remove it in your plots if you do not want it.) I do not see the intersection even when setting a wider 'LineWidth' ffor the plot3 line. The minimum of the matrix intersections is , not zero.

The surfaces might intersect if ‘Y’ (here ‘C’) is allowed to become negative for both surfaces. They do not appear to intersect with the present variable definitions.

I will help further as necessary.

.

[M]

Many thanks, @Star Strider, it is great that you plotted a plane equal to zero because in fact you kill two birds with one stone, as I wanted to know the equilibrium points of this system which is cdot=ctdot=0...Anyway, could you please tell me what is the value of ct and c that cdot(Ct,C)-ctdot(Ct,C)=1.8648e-05 which is min value?

[Star Strider]

My pleasure!

I will do my best!

I edited my previous Comment to specify:

Mtx = cdot(Ct,C)-ctdot(Ct,C);

[minMtx,idx] = min(Mtx(:))

[row,col] = ind2sub(size(Mtx), idx)

The minimum corresponds to row 1 and column 26 of ‘Mtx’ which is the matrix created by subtracting ‘ctdot’ from ‘cdot’. If you want to mark it, use either ‘idx’ or the subscripts. To determine if there is more than one value for the minimum (since min only returns the first instance it doscovers), use the find function:

idx = find(Mtx == min(Mtx(:)))

That will return all of them if there are more than one. Use the same ind2sub call I used to return the row and column values, although that would not be necessary in order to plot them. It would just tell you where they are in the matrix in a way that would be easier to interpret.

I went back and added the find call as well.

[M]

@Star Strider, I again have question about this problem but it is not related to plotting in 3d and surface, this time I wanted to plot in 2d (c-ct plane). My question is that I have to ask it here or as a new question?

[Star Strider]

It is relevant here. The simplest option would be to use the view function to rotate the plot to look at it on edge, for example:

view(90,0) % Z(Y)

view(0,0) % Z(X)

depending on what you want. Those can be used with any of the existing surf plots, or others you may want to create with these functions. (Another option is simply to use plot, however that may not be meaningful in this context.)

I have added both of the view options here —

vplc=0.16;
delta=0.1;
Ktau=0.045;
Kc=0.1; 
K=0.0075;
Kp=0.15;
gamma=5.5;
Kb=0.4;
alpha0=delta*6.81e-6/(0.002);
alpha1=delta*2.27e-5/(0.002);
Ke=7;
Vs=0.002;
Ks=0.1;
Kf=0.18;
kplc=0.055;
ki=2;
tau_max=0.176;
Vss=0.044;

ct = 0:.1:2.5;
c = 0:.1:2.5;

A=(-(Vss.*c.^2)./(Ks.^2))+((Vs.*K.*gamma.^2.*ct.^2)./(Ks.^2))+alpha0+alpha1.*((Ke.^4)./(Ke.^4+(gamma.*ct).^4));
p=(vplc.*c.^2./(c.^2+kplc.^2))./ki;
h=(-(0.4.*A.*((Kc.^4).*(Kp.^2))./((p.^2.*c.^2.*gamma.*ct.*Kf))));
G1=alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4));
G2=((1-h)./tau_max).*c.^4;
Fc=(4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2);
Fct=((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2);
Fh=(gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4);

cdot=@(ct,c)(((gamma.*Kf.*(c.^4).*(p.^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2)).*(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4)))+((gamma.*Kf.*c.^4.*p.^2.*ct)./(Kb.*Kp.^2.*Ktau.^4)).*(((1-h)./tau_max).*c.^4);
ctdot=@(ct,c)(-(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4))).*((4.*gamma.*Kf).*((c.^3.*p.^2.*h.*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2)));

[Ct,C] = meshgrid(ct,c);

figure
surf(Ct,C,cdot(Ct,C))
hold on
surf(Ct,C,zeros(size(C)), 'FaceColor','r', 'FaceAlpha',0.5)
hold off
xlabel('Ct')
ylabel('C')
title('cdot')

figure
surf(Ct,C,ctdot(Ct,C))
hold on
surf(Ct,C,zeros(size(C)), 'FaceColor','r', 'FaceAlpha',0.5)
hold off
xlabel('Ct')
ylabel('C')
title('ctdot')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
xlabel('Ct')
ylabel('C')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
view(90,0)
xlabel('Ct')
ylabel('C')
zlabel('Z')
title('Z(C)')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
view(0,0)
xlabel('Ct')
ylabel('C')
zlabel('Z')
title('Z(Ct)')

figure
surf(Ct,C,cdot(Ct,C), cat(3,ones(size(Ct)),ones(size(C)),zeros(size(Ct))), 'DisplayName','cdot')
hold on
surf(C,Ct,ctdot(Ct,C), cat(3,ones(size(Ct)),zeros(size(Ct)),ones(size(C))), 'DisplayName','ctdot')  
hold off
view(0,90)
xlabel('Ct')
ylabel('C')
zlabel('Z')
title('Surface C(Ct)')
legend('Location','best')
legend('Location','best')

figure
contour(Ct,C,cdot(Ct,C),20,'m', 'DisplayName','cdot')
hold on
contour(C,Ct,ctdot(Ct,C),20,'c', 'DisplayName','ctdot') 
hold off
% view(0,90)
xlabel('Ct')
ylabel('C')
zlabel('Z')
title('Contour  C(Ct)')
legend('Location','best')

figure
surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))
hold on
contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',5)  
hold off
xlabel('Ct')
ylabel('C')
ylim([0 0.2])
zlim([-1 1])
title('Zoomed')

Mtx = cdot(Ct,C)-ctdot(Ct,C);
[minMtx,idx] = min(Mtx(:))

minMtx = 1.8648e-05

idx = 651

[row,col] = ind2sub(size(Mtx), idx)

row = 1

col = 26

idx = find(Mtx == min(Mtx(:)))

idx = 651

[row,col] = ind2sub(size(Mtx), idx)

row = 1

col = 26

.

[M]

Thanks for the answer. The one that I want is not projection on plane, as you know ctdot and cdot are functions of c and ct; I want to plot it in a c-ct plane which might be becomes a function of c or ct, which you said "Plot" would be meaningless in this case. I attached a figure to show clearly what I ment.

[Torsten]

vplc=0.16;
delta=0.1;
Ktau=0.045;
Kc=0.1; 
K=0.0075;
Kp=0.15;
gamma=5.5;
Kb=0.4;
alpha0=delta*6.81e-6/(0.002);
alpha1=delta*2.27e-5/(0.002);
Ke=7;
Vs=0.002;
Ks=0.1;
Kf=0.18;
kplc=0.055;
ki=2;
tau_max=0.176;
Vss=0.044;
A=@(ct,c) (-(Vss.*c.^2)./(Ks.^2))+((Vs.*K.*gamma.^2.*ct.^2)./(Ks.^2))+alpha0+alpha1.*((Ke.^4)./(Ke.^4+(gamma.*ct).^4));
p=@(ct,c)(vplc.*c.^2./(c.^2+kplc.^2))./ki;
h=@(ct,c)(-(0.4.*A(ct,c).*((Kc.^4).*(Kp.^2))./((p(ct,c).^2.*c.^2.*gamma.*ct.*Kf))));
G1=@(ct,c)alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4));
G2=@(ct,c)((1-h)./tau_max).*c.^4;
Fc=@(ct,c)(4.*gamma.*Kf).*((c.^3.*p(ct,c).^2.*h(ct,c).*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2);
Fct=@(ct,c)((gamma.*Kf.*(c.^4).*(p(ct,c).^2).*h)./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2);
Fh=@(ct,c)(gamma.*Kf.*c.^4.*p(ct,c).^2.*ct)./(Kb.*Kp.^2.*Ktau.^4);
cdot=@(ct,c)(((gamma.*Kf.*(c.^4).*(p(ct,c).^2).*h(ct,c))./(Kb.*Kp.^2.*Ktau.^4))+((Vs.*K.*gamma.^2)./(Ks.^2))-((4.*gamma.^4.*ct.^3.*alpha1.*Ke.^4)./(Ke.^4+(gamma.*ct).^4).^2)).*(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4)))+((gamma.*Kf.*c.^4.*p(ct,c).^2.*ct)./(Kb.*Kp.^2.*Ktau.^4)).*(((1-h(ct,c))./tau_max).*c.^4);
ctdot=@(ct,c)(-(alpha0+(alpha1.*Ke.^4./((gamma.*ct).^4+Ke.^4))).*((4.*gamma.*Kf).*((c.^3.*p(ct,c).^2.*h(ct,c).*ct)./(Kb.*Kp.^2.*Ktau.^4))-(2.*Vss.*c./Ks.^2)));
figure(1)
fimplicit(@(x,y)cdot(x,y),[-2.5 2.5 -2.5 2.5])

figure(2)
fimplicit(@(x,y)ctdot(x,y),[-2.5 2.5 -2.5 2.5])

[Star Strider]

In the latest plots, ‘X’ is ‘Ct’ and ‘Y’ is ‘C’. Drawing them like that (essentially eliminating the ‘Z’ axis) would be looking at them from the top down, with the view(0,90) and they would appear as a flat plane. The only realistic way to plot them with ‘Ct’ as function of ‘C’ is to use a contour plot, so I added that. I doubt that anything else would work.

I edited my previous Comment to add these plots.

[M]

Thanks a million for really nice and complete expalanation. I have a question:

I labeled the contour to see the zero contour, why the curve is not the same as the last figure that you plotted?

I mean why the projection of zero-contour in xy plane for z=f(x,y) is not equal y=f(x)?

[Star Strider]

As always, my pleasure!

I am a bit lost. I am not certain what you are referring to. (The surfaces do not intersect, at least with the current definitions of ‘c’ and ‘ct’.) I plotted the contours of both functions together, since I could not see how to produce the sort of plot you described. (I still do not.)

I did not see Torsten’s comment before. Those plots are where the functions equal zero (the definition of implicit functions), however they are plotted over positive and negative values of ‘c’ and ‘ct’, neither of which are defined as being negative originally. I am not certain what to make of them in the context of this problem.

[Torsten]

I am not certain what to make of them in the context of this problem.

According to the photo @M included, these are the plots of c = g1(ct) and c = g2(ct) that are requested.

It should be no problem to restrict c and ct to positive ranges if desired.

[Star Strider]

I initially plotted them as overlapping contours of the matrix functions, since it was not obvious to me how to plot them otherwise.

This would be straightforward if the surfaces intersect, that being the original problem. They don’t, at least in the variable ranges defined.

[Torsten]

At least cdot(ct,c) = 0 is empty for ct >=0.

ctdot(ct,c) = 0 gives two implicit curves symmetric to the ct-axis.

[Star Strider]

Originally, the objective was to determine the intersection:

figure

surf(Ct,C,cdot(Ct,C)-ctdot(Ct,C))

hold on

contour3(C,Ct,cdot(Ct,C)-ctdot(Ct,C), [0 0], '-r', 'LineWidth',2)

I am not certain where we are at this point.