#!/usr/bin/env octave

disp(sprintf("#%11s %12s %12s %12s",'y_m','z_m','h_l','phi_deg'))

function sig(phi_deg,h_l)
  l_m = 32.92;
  h_m = h_l*l_m;
  y_m =  h_m*sin(phi_deg/180*pi);
  y_m = max([y_m 0.0001]);
  z_m = -h_m*cos(phi_deg/180*pi);

  disp(sprintf("%12.6f %12.6f %12.6f %12.6f",y_m,z_m,h_l,phi_deg))
end

hls=1:2:5;
phis=30:30:150;

for h_l = hls
  sig(0,h_l);
end
for phi_deg = phis
  sig(phi_deg,3);
end

disp(sprintf('VARIABLES = "X" "Y" "Z" "dp/pinf"'));

function zone(phi_deg,h_l)
  l_m = 32.92;
  h_m = h_l*l_m;
  y_m =  h_m*sin(phi_deg/180*pi);
  y_m = max([y_m 0.0001]);
  z_m = -h_m*cos(phi_deg/180*pi);

  disp(sprintf('ZONE T="HL=%3.1f,PHI=%04.1f"',h_l,phi_deg))
  disp(sprintf(" %3.1f %12.6f %12.6f %3.1f",0,y_m,z_m,0))
end

for h_l = hls
  zone(0,h_l);
end
for phi_deg = phis
  zone(phi_deg,3);
end

function nml(phi_deg,h_l,i)
  l_m = 32.92;
  h_m = h_l*l_m;
  y_m =  h_m*sin(phi_deg/180*pi);
  y_m = max([y_m 0.0001]);
  z_m = -h_m*cos(phi_deg/180*pi);

  disp(sprintf("y_ray(%2d) = %12.6f",i,y_m))
  disp(sprintf("z_ray(%2d) = %12.6f",i,z_m))
end

disp(sprintf('&sonic_boom'));
i = 0;
for h_l = hls
  i+=1;
  nml(0,h_l,i);
end
for phi_deg = phis
  i+=1;
  nml(phi_deg,3,i);
end

disp(sprintf('nsignals = %d',i));
disp(sprintf('/'));
