Parametric Plot with two parameter and a function

Question

I want to plot this parametric equation by using these equations in Matlab

p_x = p_0*coshu*cosv, p_y = p_0*sinhu*sinv

sinv*( sqrt(1 − γ)*coshu + cosα) = −sinα *sinhu

and i need a plot between p_y/p_0 vs p_x/p_0. as shown in figure

where u is a free parameter when' α = 8*pi/5. and γ = 0, 0.05, 0.15, 0.2

I tried a code solving above equations as;

close all
clear all
clc
a = 8*pi/5        % 'a' as alpha
 %z=0;           % 'z'   as  gamma
z=0.15      
u = -5:0.003:5;
x = cosh(u).*sqrt(1 - (sin(a).*sin(a).*sinh(u).*sinh(u))./square((sqrt(1-z).*cosh(u) + cos(a))));
% where x =  p_x/p_0  and y = p_y/p_0
y= -1.*((sinh(u).*sinh(u).*sin(a))./(1*sqrt(1-z).*cosh(u) + cos(a)));
plot(x,y,'-k')

Another try in Solving Equation 1 in comment(There is still somethng wrong with sign(cos(v))):

clc; clear;
alpha=8*pi/5; gamma=0.05;
t=1;
Py0={};
for Px0=-3:.5:3
  syms u
  F=cosh(u)*sqrt(1 - ((sin(alpha)*sinh(u))/(sqrt(1-gamma)*cosh(u)+cos(alpha)))^2 )-Px0;
  u=double(solve(F));
  Py0{t}=sinh(u).*(-(sin(alpha).*sinh(u))./(sqrt(1-gamma).*cosh(u)+cos(alpha)));
  t=t+1;
  clear u
end;
Py0
% plot(-3:0.5:3,Py0)

Answer 1

u is the free parameter, but its range is limited by the third equation:

sin(v)*( sqrt(1 − γ)*cosh(u) + cos(α)) = −sin(α)*sinh(u)

This can be rewritten as:

sin(v) = -sin(alpha)*sinh(u)/(sqrt(1-y)*cosh(u)+cos(alpha))

Knowning that

abs(sin(v)) <= 1

gives a condition for u.

Using that

cosh(x)^2 - sinh(x)^2 = 1

the condition becomes:

abs(-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha)) <= 1

Since cosh(x) is an even funtion, so is the expression above. Therefore it suffices to calculate

-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha) <= 1

We want to know the maximum u for which the expression hold (u_max), because then we know that u is limited in the range [-u_max,u_max]. So we need to solve

-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha) = 1

This is a second order polynomial, and will therefore have 2 solutions. We are interested in the real solutions, and if all solutions are imaginary, then there is no limit on the range of u.

Putting this in MATLAB, results in the following code:

g = [0 .05 .15 .2];         % different gammas
p = {'--k' ':g' '-r' '-b'}; % for plotting
alpha = 8*pi/5;
syms x
for i=1:length(g)
    gamma = g(i);

    % solve condition
    sinv = -sin(alpha).*sinh(x)./(sqrt(1-gamma).*cosh(x)+cos(alpha));
    sols = solve((sinv) == 1, x); % will have max 2 solutions

    % pick right solution
    if isreal(sols(1))
        u_max = double(sols(1));
    elseif isreal(sols(2))
        u_max = double(sols(2));
    else                         % both sols imaginary: no limit on u_max
        u_max = 5;
    end

    u = -u_max:0.003:u_max;
    sinv = -sin(alpha).*sinh(u)./(sqrt(1-gamma).*cosh(u)+cos(alpha));
    cosv = sqrt(1-sinv.^2); % actually +-sqrt(), taken into account when plotting

    px = cosh(u).*cosv;
    py = sinh(u).*sinv;

    plot(px,py, p{i},-px,py, p{i})
    hold on
end
hold off

---

EDIT: update of code

g = [0 .01 .1 .75];
p = {'--k' ':g' '-r' 'ob'};
alpha = 4*pi/3;
syms x
for i=1:4
    gamma = g(i);
    interval = inf;
    sinv = -sin(alpha).*sinh(x)./(sqrt(1-gamma).*cosh(x)+cos(alpha));
    sols = solve((sinv) == 1, x); % will have max 2 solutions

    if length(sols) > 1
        if isreal(sols(1)) && isreal(sols(2))  % if there are 2 real solutions, interval between is valid or unvalid
            if eval(subs(sinv,x,(double(sols(1))+double(sols(2)))/2)) > 1 %interval inbetween is unvalid => u ok everywhere except in interval
                u_max = double(min(sols(1),sols(2)));
                u_min = double(max(sols(1),sols(2)));
                interval = 0;
            else  %interval inbetween is valid => u ok in interval
                u_max = double(max(sols(1),sols(2)));
                u_min = double(min(sols(1),sols(2)));
                interval = 1;
            end
        elseif isreal(sols(1))
            u_max = double(sols(1));
        elseif isreal(sols(2))
            u_max = double(sols(2));
        else
            u_max = 3;
        end
    elseif isreal(sols)
        if eval(subs(sinv,x,sols-.1)) < 1 && eval(subs(sinv,x,sols+.1)) < 1
            u_max = 3;
        else
            u_max = double(sols);
        end
    elseif eval(subs(sinv,x,1)) < 1
        u_max = 3;
    else
        u_max = 0;
    end

    if interval == 1
        u1 = u_min:0.003:u_max;
        u2 = -u1;
        sinv1 = -sin(alpha).*sinh(u1)./(sqrt(1-gamma).*cosh(u1)+cos(alpha));
        sinv2 = -sin(alpha).*sinh(u2)./(sqrt(1-gamma).*cosh(u2)+cos(alpha));
        cosv1 = sqrt(1-sinv1.^2);
        cosv2 = sqrt(1-sinv2.^2);


        if imag(cosv1) < 10^(-6)
            cosv1 = real(cosv1);
        end
        if imag(cosv2) < 10^(-6)
            cosv2 = real(cosv2);
        end
        if imag(cosv1) < 10^(-6)
            cosv1 = real(cosv1);
        end
        if imag(cosv2) < 10^(-6)
            cosv2 = real(cosv2);
        end

        px1 = cosh(u1).*cosv1;
        py1 = sinh(u1).*sinv1;
        px2 = cosh(u2).*cosv2;
        py2 = sinh(u2).*sinv2;

        plot(([-px1(end:-1:1) px1]),([py1(end:-1:1) py1]), p{i}, ([-px2(end:-1:1) px2]),([py2(end:-1:1) py2]), p{i})
        hold on
    elseif interval == 0
        u1 = -u_max:0.003:u_max;
        u2 = u_min:0.003:3;
        sinv1 = -sin(alpha).*sinh(u1)./(sqrt(1-gamma).*cosh(u1)+cos(alpha));
        sinv2 = -sin(alpha).*sinh(u2)./(sqrt(1-gamma).*cosh(u2)+cos(alpha));
        cosv1 = sqrt(1-sinv1.^2);
        cosv2 = sqrt(1-sinv2.^2);

        if imag(cosv1) < 10^(-6)
            cosv1 = real(cosv1);
        end
        if imag(cosv2) < 10^(-6);
            cosv2 = real(cosv2);
        end

        px1 = cosh(u1).*cosv1;
        py1 = sinh(u1).*sinv1;
        px2 = cosh(u2).*cosv2;
        py2 = sinh(u2).*sinv2;

        plot([-px1(end:-1:1) (px1)],[py1(end:-1:1) (py1)], p{i}, ([-px2(end:-1:1) px2]),([py2(end:-1:1) py2]), p{i})
        hold on
    else
        u = -u_max:0.003:u_max;
        sinv1 = -sin(alpha).*sinh(u)./(sqrt(1-gamma).*cosh(u)+cos(alpha));
        cosv1 = sqrt(1-sinv1.^2);

         if imag(cosv1) < 10^(-6)
            cosv1 = real(cosv1);
        end

        px1 = cosh(u).*cosv1;
        py1 = sinh(u).*sinv1;

        plot(-px1(end:-1:1), py1(end:-1:1), p{i}, px1, py1, p{i})
        hold on
    end 
end
hold off