Trigonometry

% ──── Angle Conversion in GNU Octave ────────────────────────────
% Octave works natively in RADIANS for trig functions

% Degrees to radians
degrees = [0, 30, 45, 60, 90, 180, 360];
radians = degrees .* (pi ./ 180);
printf("Degrees: "); disp(degrees);
printf("Radians: "); disp(radians);

% Octave built-in conversion
r1 = deg2rad(45);     % 0.7854
d1 = rad2deg(pi/4);   % 45.000

% Arc length and sector area
r = 5;           % radius = 5 units
theta = pi/3;   % 60 degrees
arc_len  = r * theta;              % s = r * θ
sec_area = 0.5 * r^2 * theta;     % A = ½r²θ
printf("Arc length : %.4f\n", arc_len);
printf("Sector area: %.4f\n", sec_area);

% Visualise angle conversion table
printf("\n%-10s %-15s %-15s\n", "Degrees", "Radians", "Fraction of π");
for d = [0,30,45,60,90,120,180,270,360]
  printf("%-10d %-15.6f %-15.6f\n", d, deg2rad(d), deg2rad(d)/pi);
endfor

% ──── Plot the Unit Circle with special angles ───────────────────
figure('Name', 'Unit Circle');
theta_fine = 0:0.01:2*pi;
plot(cos(theta_fine), sin(theta_fine), 'b-', 'LineWidth', 2); hold on;

% Draw axes
line([-1.3 1.3], [0 0], 'Color', 'k');
line([0 0], [-1.3 1.3], 'Color', 'k');

% Plot special angles
angles = [0, pi/6, pi/4, pi/3, pi/2, 2*pi/3, 3*pi/4, 5*pi/6, pi, ...
           7*pi/6, 5*pi/4, 4*pi/3, 3*pi/2, 5*pi/3, 7*pi/4, 11*pi/6];
labels = {'0','π/6','π/4','π/3','π/2','2π/3','3π/4','5π/6','π', ...
          '7π/6','5π/4','4π/3','3π/2','5π/3','7π/4','11π/6'};

for i = 1:length(angles)
  cx = cos(angles(i)); cy = sin(angles(i));
  plot(cx, cy, 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
  line([0 cx], [0 cy], 'Color', [0.7 0.7 0.7], 'LineStyle', '--');
  text(cx*1.15, cy*1.15, labels{i}, 'FontSize', 9, 'HorizontalAlignment', 'center');
endfor

axis equal; axis([-1.5 1.5 -1.5 1.5]);
title('The Unit Circle'); xlabel('cos θ'); ylabel('sin θ');
grid on; hold off;

% ──── Six Trig Functions in GNU Octave ───────────────────────────
theta_deg = 30;
theta     = deg2rad(theta_deg);

s  = sin(theta);        % sine
c  = cos(theta);        % cosine
t  = tan(theta);        % tangent
cs = 1/sin(theta);     % cosecant
se = 1/cos(theta);     % secant
co = 1/tan(theta);     % cotangent

printf("θ = %d° = %.4f rad\n", theta_deg, theta);
printf("sin = %.6f\ncos = %.6f\ntan = %.6f\n", s, c, t);
printf("csc = %.6f\nsec = %.6f\ncot = %.6f\n", cs, se, co);

% ──── Table for all key angles ────────────────────────────────────
printf("\n%-6s %-10s %-10s %-10s\n", "Deg", "sin", "cos", "tan");
for d = [0,30,45,60,90,120,135,150,180]
  r = deg2rad(d);
  if abs(cos(r)) < 1e-10
    printf("%-6d %-10.6f %-10.6f %-10s\n", d, sin(r), cos(r), "undef");
  else
    printf("%-6d %-10.6f %-10.6f %-10.6f\n", d, sin(r), cos(r), tan(r));
  endif
endfor

% ──── Verify Pythagorean identity: sin²+cos² = 1 ─────────────────
for d = [0:15:360]
  r   = deg2rad(d);
  val = sin(r)^2 + cos(r)^2;
  assert(val, 1, 1e-12);   % should not throw
endfor
disp("Pythagorean identity verified for all multiples of 15°");

% ──── Right Triangle Solver ───────────────────────────────────────
function solve_right_triangle(theta_deg, hyp)
  theta = deg2rad(theta_deg);
  opp   = hyp * sin(theta);
  adj   = hyp * cos(theta);
  printf("θ=%.1f°  hyp=%.4f  opp=%.4f  adj=%.4f\n", theta_deg, hyp, opp, adj);
endfunction

solve_right_triangle(35, 10);
solve_right_triangle(60, 5);

% ──── Law of Cosines: find side c given a, b, C ───────────────────
function c = law_of_cosines(a, b, C_deg)
  C = deg2rad(C_deg);
  c = sqrt(a^2 + b^2 - 2*a*b*cos(C));
endfunction

c = law_of_cosines(7, 10, 45);
printf("Law of Cosines: c = %.6f\n", c);

% ──── Law of Sines: find angle B given a, b, A ────────────────────
function B_deg = law_of_sines_angle(a, b, A_deg)
  A   = deg2rad(A_deg);
  B   = asin(b * sin(A) / a);
  B_deg = rad2deg(B);
endfunction

B = law_of_sines_angle(8, 5, 50);
printf("Law of Sines: B = %.4f°\n", B);

% ──── Heron's Formula ────────────────────────────────────────────
function A = herons_area(a, b, c)
  s = (a + b + c) / 2;
  A = sqrt(s * (s-a) * (s-b) * (s-c));
endfunction

area = herons_area(5, 7, 9);
printf("Heron area of (5,7,9) triangle: %.6f\n", area);

% ──── Two-sided area formula ──────────────────────────────────────
area2 = 0.5 * 5 * 7 * sin(deg2rad(60));
printf("½ab·sinC area (a=5,b=7,C=60°): %.6f\n", area2);

% ──── Verify trig identities numerically ─────────────────────────
t = pi/7;   % arbitrary test angle
A = pi/5;
B = pi/8;
tol = 1e-12;

% Pythagorean
assert(sin(t)^2 + cos(t)^2,     1,             tol);
assert(tan(t)^2 + 1,              1/cos(t)^2,    tol);
disp("Pythagorean identities: OK");

% Double angle
assert(sin(2*t), 2*sin(t)*cos(t), tol);
assert(cos(2*t), cos(t)^2 - sin(t)^2, tol);
disp("Double angle identities: OK");

% Sum identities
assert(sin(A+B), sin(A)*cos(B) + cos(A)*sin(B), tol);
assert(cos(A+B), cos(A)*cos(B) - sin(A)*sin(B), tol);
disp("Sum/difference identities: OK");

% Product-to-sum
assert(sin(A)*sin(B), 0.5*(cos(A-B) - cos(A+B)), tol);
assert(cos(A)*cos(B), 0.5*(cos(A-B) + cos(A+B)), tol);
disp("Product-to-sum identities: OK");

% Half-angle (for positive quadrant only)
t2 = pi/5;
assert(sin(t2/2), sqrt((1 - cos(t2))/2), tol);
disp("Half-angle identity: OK");

% ──── Inverse Trig Functions in Octave ───────────────────────────
% asin, acos, atan return radians
x = 0.5;
a1 = asin(x);          % = π/6
a2 = acos(x);          % = π/3
a3 = atan(x);          % ≈ 0.4636

printf("asin(0.5) = %.6f rad = %.4f°\n", a1, rad2deg(a1));
printf("acos(0.5) = %.6f rad = %.4f°\n", a2, rad2deg(a2));
printf("atan(0.5) = %.6f rad = %.4f°\n", a3, rad2deg(a3));

% atan2 — four-quadrant version
printf("\natan2 examples:\n");
printf("atan2(1, 1)   = %.4f° (Q1)\n", rad2deg(atan2(1,1)));
printf("atan2(1,-1)   = %.4f° (Q2)\n", rad2deg(atan2(1,-1)));
printf("atan2(-1,-1)  = %.4f° (Q3)\n", rad2deg(atan2(-1,-1)));
printf("atan2(-1, 1)  = %.4f° (Q4)\n", rad2deg(atan2(-1,1)));

% Verify: arcsin(x) + arccos(x) = π/2
for x = -1:0.25:1
  assert(asin(x) + acos(x), pi/2, 1e-12);
endfor
disp("arcsin(x) + arccos(x) = π/2 verified");

% Plot inverse functions
figure('Name', 'Inverse Trig');
x1 = linspace(-1, 1, 500);
x2 = linspace(-5, 5, 500);
subplot(1,3,1); plot(x1, rad2deg(asin(x1))); title('arcsin'); grid on;
subplot(1,3,2); plot(x1, rad2deg(acos(x1))); title('arccos'); grid on;
subplot(1,3,3); plot(x2, rad2deg(atan(x2))); title('arctan'); grid on;

% ──── Sinusoidal wave transformations ────────────────────────────
x = linspace(0, 4*pi, 1000);

% y = A·sin(Bx + C) + D
A = 3;  B = 2;  C = pi/4;  D = 1;
y1 = A * sin(B*x + C) + D;
y2 = sin(x);                  % reference
y3 = cos(x);

figure('Name', 'Wave Transformations');
plot(x, y2, 'b-', 'LineWidth', 1.5); hold on;
plot(x, y3, 'r-', 'LineWidth', 1.5);
plot(x, y1, 'g-', 'LineWidth', 2);
legend('sin(x)', 'cos(x)', '3sin(2x+π/4)+1');
title('Sinusoidal Transformations'); grid on; hold off;

% ──── Sound wave: 440 Hz A note ──────────────────────────────────
Fs = 44100;          % sample rate Hz
dur = 1;             % 1 second
t  = 0:1/Fs:dur;
f  = 440;             % 440 Hz = concert A
wave = sin(2*pi*f*t);

figure('Name', '440 Hz A Note');
plot(t(1:500), wave(1:500));  % show first ~11ms
title('440 Hz sine wave (Concert A)');
xlabel('Time (s)'); ylabel('Amplitude'); grid on;

% ──── Phase shift visualisation ──────────────────────────────────
x = linspace(0, 2*pi, 500);
figure('Name', 'Phase Shifts');
phases = [0, pi/6, pi/4, pi/3, pi/2];
colours = {'b', 'r', 'g', 'm', 'c'};
hold on;
for i = 1:length(phases)
  plot(x, sin(x + phases(i)), colours{i}, 'LineWidth', 1.5);
endfor
legend('φ=0','φ=π/6','φ=π/4','φ=π/3','φ=π/2');
title('Phase Shifted Sine Waves'); grid on; hold off;

% ──── Polar Curves ────────────────────────────────────────────────
theta = linspace(0, 2*pi, 1000);

% Cardioid: r = 1 + cos(θ)
r_card = 1 + cos(theta);
figure('Name', 'Polar Curves');
subplot(2,3,1);
polar(theta, r_card); title('Cardioid: r=1+cos θ');

% Rose curve (4 petals): r = cos(2θ)
r_rose = cos(2*theta);
subplot(2,3,2);
polar(theta, r_rose); title('Rose: r=cos(2θ)');

% Archimedean spiral: r = θ/(2π)
theta_sp = linspace(0, 6*pi, 1000);
r_sp = theta_sp / (2*pi);
subplot(2,3,3);
polar(theta_sp, r_sp); title('Spiral: r=θ/2π');

% ──── Euler's formula: visualise e^(iθ) path ─────────────────────
theta_e = linspace(0, 2*pi, 1000);
z = exp(1i * theta_e);         % complex unit circle
subplot(2,3,4);
plot(real(z), imag(z)); axis equal; grid on;
title('e^{iθ} in complex plane');

% Verify Euler's formula
t = pi/4;
lhs = exp(1i * t);
rhs = cos(t) + 1i*sin(t);
assert(abs(lhs - rhs), 0, 1e-14);
printf("Euler's formula verified: e^(iπ/4) = %.6f + %.6fi\n", real(lhs), imag(lhs));

% De Moivre's theorem: (cos θ + i sin θ)^n = cos(nθ) + i sin(nθ)
n = 5; t = pi/7;
lhs = (cos(t) + 1i*sin(t))^n;
rhs = cos(n*t) + 1i*sin(n*t);
printf("De Moivre error: %.2e\n", abs(lhs - rhs));

% Polar to Cartesian and back
r0 = 5; theta0 = pi/3;
x0 = r0 * cos(theta0);
y0 = r0 * sin(theta0);
r_check = sqrt(x0^2 + y0^2);
t_check = atan2(y0, x0);
printf("Round-trip: r=%.4f (orig %.4f), θ=%.4f (orig %.4f)\n", ...
       r_check, r0, t_check, theta0);

% ──── Fourier Series: Square Wave Approximation ───────────────────
% A square wave can be built from odd harmonics of sine
x = linspace(0, 2*pi, 2000);
N_terms = [1, 3, 7, 15, 51];
figure('Name', 'Fourier Series');
for idx = 1:length(N_terms)
  N = N_terms(idx);
  y = zeros(1, length(x));
  for k = 1:2:N       % odd harmonics only
    y = y + (4/(pi*k)) * sin(k*x);
  endfor
  subplot(length(N_terms), 1, idx);
  plot(x, y); ylim([-1.5, 1.5]);
  title(sprintf('Square wave: %d terms', ceil(N/2)));
endfor

% ──── FFT Analysis ────────────────────────────────────────────────
Fs = 1000;                 % sampling frequency
T  = 1/Fs;                 % sample period
L  = 1000;                % signal length
t  = (0:L-1)*T;

% Signal = 50 Hz + 120 Hz components
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

Y  = fft(S);              % Octave FFT
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f  = Fs*(0:L/2)/L;

figure('Name', 'FFT');
subplot(2,1,1); plot(t(1:200), S(1:200));
title('Signal: 50 Hz + 120 Hz'); xlabel('t (s)'); grid on;
subplot(2,1,2); plot(f, P1);
title('Single-Sided Spectrum'); xlabel('f (Hz)'); ylabel('|P1(f)|'); grid on;

% ──── Discrete Cosine Transform (DCT) — basis of JPEG ────────────
N = 8;
function C = dct_matrix(N)
  C = zeros(N, N);
  for k = 0:N-1
    for n = 0:N-1
      if k == 0
        C(k+1,n+1) = 1/sqrt(N);
      else
        C(k+1,n+1) = sqrt(2/N) * cos(pi*k*(2*n+1)/(2*N));
      endif
    endfor
  endfor
endfunction

D = dct_matrix(8);
x_block = [52 55 61 66 70 61 64 73]; % 8-pixel block
X_dct   = D * x_block';
printf("DCT of 8-pixel block:\n"); disp(X_dct');

% ──── Linear Congruential Generator ─────────────────────────────
function seq = lcg(seed, a, c, m, N)
  seq = zeros(1, N);
  x   = seed;
  for i = 1:N
    x      = mod(a*x + c, m);
    seq(i) = x;
  endfor
endfunction

% Numerical Recipes parameters (32-bit)
a = 1664525; c = 1013904223; m = 2^32;
seq = lcg(42, a, c, m, 20);
normalised = seq / m;       % uniform [0,1)
printf("LCG sequence (first 5, normalised): ");
disp(normalised(1:5));

% ──── Sinusoidal XOR Stream Cipher (educational demo) ────────────
% NOTE: This is NOT a secure cipher — for illustration only!

function cipher_text = sin_xor_encrypt(plaintext_bytes, key_freq, key_phase)
  n           = length(plaintext_bytes);
  t           = (0:n-1);
  keystream    = uint8(mod(floor(128 * (1 + sin(key_freq*t + key_phase))), 256));
  cipher_text  = bitxor(plaintext_bytes, keystream);
endfunction

msg     = uint8('Hello, Trig!');
freq    = 0.37;    % key frequency
phase   = 1.23;    % key phase

ciphertext  = sin_xor_encrypt(msg, freq, phase);
decrypted   = sin_xor_encrypt(ciphertext, freq, phase);  % XOR is its own inverse

printf("Original  : %s\n", char(msg));
printf("Encrypted : "); disp(ciphertext);
printf("Decrypted : %s\n", char(decrypted));
assert(all(decrypted == msg));
disp("Decryption verified!");

% ──── Diffie-Hellman Key Exchange Simulation ─────────────────────
function result = powmod(base, exp_val, modulus)
  % Modular exponentiation: base^exp_val mod modulus
  result = 1;
  base   = mod(base, modulus);
  while exp_val > 0
    if mod(exp_val, 2) == 1
      result = mod(result * base, modulus);
    endif
    exp_val = floor(exp_val / 2);
    base    = mod(base * base, modulus);
  endwhile
endfunction

% Public parameters
p = 23;   % prime modulus (small for demo; real DH uses 2048+ bit primes)
g = 5;    % generator (primitive root mod p)

% Alice's secret key
a = 6;
A = powmod(g, a, p);   % g^a mod p  (Alice sends this publicly)

% Bob's secret key
b = 15;
B = powmod(g, b, p);   % g^b mod p  (Bob sends this publicly)

% Compute shared secret
K_alice = powmod(B, a, p);   % B^a mod p
K_bob   = powmod(A, b, p);   % A^b mod p

printf("p=%d, g=%d\n", p, g);
printf("Alice: secret a=%d, public A=%d\n", a, A);
printf("Bob  : secret b=%d, public B=%d\n", b, B);
printf("Shared secret: Alice=%d, Bob=%d\n", K_alice, K_bob);
assert(K_alice, K_bob);
disp("Shared secrets match!");

% ──── RSA Encryption / Decryption (educational toy key sizes) ────
function result = powmod(base, exp_val, modulus)
  result = 1; base = mod(base, modulus);
  while exp_val > 0
    if mod(exp_val,2)==1; result=mod(result*base,modulus); endif
    exp_val=floor(exp_val/2); base=mod(base*base,modulus);
  endwhile
endfunction

function inv = modinv(a, m)
  % Extended Euclidean algorithm
  [g, x, ~] = gcd(a, m);
  if g != 1; error('No inverse'); endif
  inv = mod(x, m);
endfunction

% Key generation
p = 61;  q = 53;     % small primes — use 1024+ bit primes in real RSA
n = p * q;                % n = 3233
phi_n = (p-1) * (q-1);  % Euler's totient = 3120
e = 17;                   % public exponent, gcd(e,φ(n)) must be 1
d = modinv(e, phi_n);     % private exponent

printf("RSA Keys:\n  n=%d  φ(n)=%d  e=%d  d=%d\n", n, phi_n, e, d);

% Encrypt / Decrypt a single character
M = 65;                   % 'A' in ASCII
C = powmod(M, e, n);      % ciphertext: M^e mod n
R = powmod(C, d, n);      % recovered:  C^d mod n

printf("Plaintext M=%d  Ciphertext C=%d  Recovered=%d\n", M, C, R);
assert(R, M);
disp("RSA encrypt/decrypt verified!");

% Verify Euler's theorem: M^φ(n) ≡ 1 (mod n) when gcd(M,n)=1
assert(powmod(M, phi_n, n), 1);
disp("Euler's theorem verified: M^φ(n) ≡ 1 (mod n)");

% ──── Elliptic Curve Point Addition (over ℝ, for visualisation) ──
function [x3, y3] = ec_add(x1, y1, x2, y2, a)
  % Add two points on y² = x³ + ax + b
  if x1 == x2
    m = (3*x1^2 + a) / (2*y1);   % tangent (point doubling)
  else
    m = (y2 - y1) / (x2 - x1);    % secant slope
  endif
  x3 = m^2 - x1 - x2;
  y3 = m*(x1 - x3) - y1;
endfunction

% Curve: y² = x³ − x (a=−1, b=0)
a = -1;  b = 0;
x_curve = linspace(-1.5, 2, 500);
figure('Name', 'Elliptic Curve');
y_sq = x_curve.^3 + a*x_curve + b;
y_pos = sqrt(max(y_sq, 0));
plot(x_curve, y_pos, 'b-', x_curve, -y_pos, 'b-');
hold on;

% Two points on the curve
x1 = -0.5; y1 = sqrt(x1^3 + a*x1 + b);
x2 = 1.5;  y2 = sqrt(x2^3 + a*x2 + b);
[x3, y3] = ec_add(x1, y1, x2, y2, a);

plot([x1 x2 x3], [y1 y2 -y3], 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
plot([x1 x3], [y1 y2], 'g--');
title('Elliptic Curve y²=x³−x: Point Addition');
grid on; axis([-1.8 2.2 -3 3]); hold off;

% ──── CORDIC Algorithm for sin/cos ───────────────────────────────
function [s, c] = cordic(angle_rad, iterations)
  % Precompute arctangent table
  atan_table = atan(2.^(-(0:iterations-1)));

  % CORDIC gain factor
  K = prod(cos(atan_table));

  x = 1; y = 0; z = angle_rad;

  for i = 0:iterations-1
    d = 1;
    if z < 0; d = -1; endif
    x_new = x - d * y * 2^(-i);
    y_new = y + d * x * 2^(-i);
    z_new = z - d * atan_table(i+1);
    x = x_new; y = y_new; z = z_new;
  endfor

  c = K * x;   % cosine
  s = K * y;   % sine
endfunction

angles = [0, pi/6, pi/4, pi/3, pi/2];
printf("\nCORDIC vs Built-in (16 iterations):\n");
printf("%-10s %-12s %-12s %-12s %-12s\n", "Angle", "CORDIC sin", "Exact sin", "CORDIC cos", "Exact cos");
for a = angles
  [s, c] = cordic(a, 16);
  printf("%-10.4f %-12.8f %-12.8f %-12.8f %-12.8f\n", a, s, sin(a), c, cos(a));
endfor

% ──── trig_toolkit.m — all six functions + inverse + checks ──────
function trig_report(theta_deg)
  t  = deg2rad(theta_deg);
  printf("─────────────────────────────────────────────\n");
  printf("θ = %g° = %.6f rad\n", theta_deg, t);
  printf("sin   = %+.8f\n", sin(t));
  printf("cos   = %+.8f\n", cos(t));
  if abs(cos(t)) > 1e-10
    printf("tan   = %+.8f\n", tan(t));
    printf("sec   = %+.8f\n", 1/cos(t));
  else
    printf("tan   = undefined\n");
    printf("sec   = undefined\n");
  endif
  if abs(sin(t)) > 1e-10
    printf("csc   = %+.8f\n", 1/sin(t));
    printf("cot   = %+.8f\n", cos(t)/sin(t));
  else
    printf("csc   = undefined\n");
    printf("cot   = undefined\n");
  endif
  printf("sin²+cos² = %.15f\n", sin(t)^2+cos(t)^2);
  printf("Quadrant  = %d\n", ceil(mod(theta_deg,360)/90 + 1e-9));
endfunction

trig_report(30);
trig_report(135);
trig_report(210);
trig_report(315);

% ──── Haversine Formula: great-circle distance on a sphere ───────
function d = haversine(lat1, lon1, lat2, lon2)
  % Inputs in degrees; output in km
  R  = 6371;            % Earth radius km
  phi1 = deg2rad(lat1); phi2 = deg2rad(lat2);
  dphi = deg2rad(lat2 - lat1);
  dlam = deg2rad(lon2 - lon1);
  a = sin(dphi/2)^2 + cos(phi1)*cos(phi2)*sin(dlam/2)^2;
  c = 2 * atan2(sqrt(a), sqrt(1-a));
  d = R * c;
endfunction

% New York (40.71°N, 74.01°W) to London (51.51°N, 0.13°W)
d = haversine(40.71, -74.01, 51.51, -0.13);
printf("NYC → London: %.1f km\n", d);   % ~5570 km

% ──── Bearing between two coordinates ────────────────────────────
function bearing = initial_bearing(lat1, lon1, lat2, lon2)
  phi1 = deg2rad(lat1); phi2 = deg2rad(lat2);
  dlam = deg2rad(lon2 - lon1);
  y    = sin(dlam) * cos(phi2);
  x    = cos(phi1)*sin(phi2) - sin(phi1)*cos(phi2)*cos(dlam);
  bearing = mod(rad2deg(atan2(y, x)), 360);
endfunction

b = initial_bearing(40.71, -74.01, 51.51, -0.13);
printf("NYC → London bearing: %.2f°\n", b);   % ~51°

% ──── Lissajous Figures: parametric x=sin(at+δ), y=sin(bt) ──────
t = linspace(0, 2*pi, 5000);
params = [1,1; 1,2; 1,3; 2,3; 3,4; 3,5];
delta = pi/4;
figure('Name', 'Lissajous Figures');
for i = 1:6
  a = params(i,1); b = params(i,2);
  x = sin(a*t + delta);
  y = sin(b*t);
  subplot(2,3,i);
  plot(x, y, 'b-', 'LineWidth', 0.8);
  axis equal; axis([-1.2 1.2 -1.2 1.2]);
  title(sprintf('a=%d, b=%d', a, b));
  grid on;
endfor

% ──── Torus: parametric surface using trig ────────────────────────
R = 3;  r = 1;      % major and minor radii
u = linspace(0, 2*pi, 80);
v = linspace(0, 2*pi, 80);
[U, V] = meshgrid(u, v);

X = (R + r*cos(V)) .* cos(U);
Y = (R + r*cos(V)) .* sin(U);
Z = r * sin(V);

figure('Name', 'Torus');
surf(X, Y, Z, 'EdgeAlpha', 0.1); colormap hsv;
title('Torus: R=3, r=1'); axis equal; shading interp;

% ──── Sphere ─────────────────────────────────────────────────────
theta_s = linspace(-pi/2, pi/2, 60);
phi_s   = linspace(-pi, pi, 60);
[TH, PH] = meshgrid(theta_s, phi_s);
Xs = cos(TH) .* cos(PH);
Ys = cos(TH) .* sin(PH);
Zs = sin(TH);
figure('Name', 'Sphere');
surf(Xs, Ys, Zs); axis equal; title('Unit Sphere'); shading interp;

% ──── Numerical solution of trig equations ───────────────────────
% Solve: sin(x) = 0.6  in [0, 2π]
val = 0.6;
x1  = asin(val);                 % principal solution in [−π/2, π/2]
x2  = pi - x1;                   % second solution in [0, 2π]
printf("sin(x)=0.6  →  x = %.4f° or %.4f°\n", rad2deg(x1), rad2deg(x2));

% Solve: 2cos²(x) − cos(x) − 1 = 0  (substitution u=cos x)
% → (2u+1)(u−1) = 0  → u = −½ or u = 1
sol1_deg = rad2deg(acos(-0.5));  % 120°
sol2_deg = rad2deg(acos(1));      %   0°
printf("2cos²x−cosx−1=0 →  x = %.1f°, %.1f° (and coterminals)\n", sol1_deg, sol2_deg);

% ──── Plot to verify ─────────────────────────────────────────────
x = linspace(0, 2*pi, 1000);
y = 2*cos(x).^2 - cos(x) - 1;
figure('Name', 'Trig Equation');
plot(x, y, 'b-', 'LineWidth', 2);
hold on; yline(0, 'r--');
title('2cos²x − cos x − 1 = 0'); xlabel('x (rad)'); grid on; hold off;

% ──── Simplified PLL: sin × sin product as phase detector ────────
% PD output: sin(θ_ref) × sin(θ_vco) = ½[cos(θ_ref−θ_vco) − cos(θ_ref+θ_vco)]
% Lowpass filter removes the sum term; ½cos(Δθ) drives VCO correction
Fs      = 1e5;
t       = 0:1/Fs:0.01;
f_ref   = 1000;
f_vco   = 1020;    % slight offset
ref_sig = sin(2*pi*f_ref*t);
vco_sig = sin(2*pi*f_vco*t);
pd_out  = ref_sig .* vco_sig;   % phase detector

% Simple moving-average LPF (approximate)
win = round(Fs / f_ref / 4);
lp  = conv(pd_out, ones(1,win)/win, 'same');

figure('Name', 'PLL Demo');
subplot(3,1,1); plot(t(1:500), ref_sig(1:500)); title('Reference 1kHz'); grid on;
subplot(3,1,2); plot(t(1:500), vco_sig(1:500)); title('VCO 1020Hz'); grid on;
subplot(3,1,3); plot(t, lp); title('Phase Detector LPF Output (∝ cos Δθ)'); grid on;

% ──── Spherical Law of Cosines ────────────────────────────────────
% For a spherical triangle with sides a,b,c (arcs) and angles A,B,C:
% cos(c) = cos(a)cos(b) + sin(a)sin(b)cos(C)
function c = spherical_law_of_cosines(a_deg, b_deg, C_deg)
  a = deg2rad(a_deg); b = deg2rad(b_deg); C = deg2rad(C_deg);
  c = acos(cos(a)*cos(b) + sin(a)*sin(b)*cos(C));
  c = rad2deg(c);
endfunction

% Triangle on Earth's surface:
% Pole → Equator/0° → Equator/90°E
c_arc = spherical_law_of_cosines(90, 90, 90);
printf("Spherical LOC: third arc = %.4f°\n", c_arc); % should be 90°

% ──── Solar Altitude Angle ────────────────────────────────────────
function alt = solar_altitude(lat_deg, decl_deg, hour_angle_deg)
  lat    = deg2rad(lat_deg);
  decl   = deg2rad(decl_deg);
  ha     = deg2rad(hour_angle_deg);
  sin_alt = sin(lat)*sin(decl) + cos(lat)*cos(decl)*cos(ha);
  alt     = rad2deg(asin(sin_alt));
endfunction

% Solar altitude at noon (hour_angle=0) at latitude 40°N on summer solstice (decl=23.5°)
alt_noon = solar_altitude(40, 23.5, 0);
printf("Solar altitude at noon (lat=40°N, summer solstice): %.2f°\n", alt_noon);