{~ Conversion of acoustic interval into vehicle velocity ~}

program Doppl2v{er tone height change to velocity, by Jenik Hollan};
uses {crt,}str_num;
const
    cl=#13#10;
    base: word = 440; del:byte=3;
    c:real=340; w:real=0; b:real=0; ht2:real=0;
    xmin:real=5;
var ht,x,a,v,vr,fr,xs:real;
    wrong:integer;
    i:word;

procedure help;
 begin
   writeln(
 'Doppl2v  <change in tone height (6 tones = 1 octave)>',cl,
 '  [ <rest frame sound frequency (default 440)/ 1 Hz>',cl,
 '  [ <least distance from you (default 5m)/ 1 m>',cl,
 '  [ <wind velocity/(1m/s)>:<azimuth/1degree from the vehicle apex> ]]]',cl,
 'says velocity of the passing vehicle.',cl2,

'( (C) Jan Hollan, N.Copernicus Observatory and Planetarium in Brno, 1999;'+cl+
' subject to the GNU General Public License, http://www.gnu.org/copyleft;'+cl+
' source code available at http://astro.sci.muni.cz/pub/hollan/programmes)'
);
   halt;
 end;

procedure showsound;
 begin
 {
  xs:=del*v/1000;
  x:=del*v;
  for i:=1 to 999 do
   begin
    x:=x-xs;
    a:=arctan(xmin/x);
    vr:=1/(1-v*cos(a)/(c+w*cos(a+b)));
    sound(round(base*vr));
    delay(del);
   end;
    sound(round(base*vr));
    delay(2*del);
  x:=del*v/1000;
  for i:=1 to 999 do
   begin
    x:=x+xs;
    a:=arctan(xmin/x);
    vr:=1/(1+v*cos(a)/(c-w*cos(b-a)));
    sound(round(base*vr));
    delay(del);
   end;
    nosound;
  }
 end;

begin
 if paramcount>0 then
  begin
   ht:=R_S(ItemStrD(':',1,paramstr(1)));
   if ht=0 then help;
   if ItemCountD(':',paramstr(1))=2 then
    ht2:=R_S(ItemStrD(':',2,paramstr(1)));
   if Paramcount>1 then
    begin
     val(paramstr(2),i,wrong);
     if wrong=0 then base:=i;
     if Paramcount>2 then
      begin
       val(paramstr(3),v,wrong);
       if wrong=0 then xmin:=v;
      end;
     if Paramcount>3 then
      begin
       if ItemCountD(':',paramstr(4))=2 then
        begin
         v:=R_S(ItemStrD(':',1,paramstr(4)));
         b:=pi*R_S(ItemStrD(':',2,paramstr(4)))/180;
         w:=v*cos(b);
         a:=(c-w)/(c+w);
        end;
      end;
    end
  end
 else help;
 fr:=exp(ht*ln(2)/6);
 if w=0 then
  v:=c*(1 - 2/(fr+1))
 else
  v:=(c-w)*(1 - (1+a)/(fr+a));
 writeln('Tone height change amounting to h= ',ht:4:1,' tone',cl,
 ' (6 tones = 12 halftones being one octave)');
 writeln('implies the vehicle passed you at ',
    v*3.6:6:1,' km/h ',cl,
         '                                = ',v:6:1,' m/s');
 if w=0 then
  writeln(
    ' using v = c*(1-2/(2^(h/6)+1))',cl,
    ' and c=340 m/s as the velocity of sound.')
 else
  writeln(
' using v= (c-w)*(1-(1+a)/(2^(h/6)+a))',cl,
' where w is velocity of the component of wind parallel to vehicle velocity',cl,
' a = (c-w) / (c+w)',cl,
' and c = 340 m/s is the velocity of sound.');

if ht2>0 then
 begin
  writeln(cl,'interval   velocity',cl,
             '--------  ----  ----',cl,
             ' 1 tone   km/h   m/s',cl);
  while ht<=ht2 do
   begin
    fr:=exp(ht*ln(2)/6);
    if w=0 then
     v:=c*(1 - 2/(fr+1))
    else
     v:=(c-w)*(1 - (1+a)/(fr+a));
    writeln(ht:5:1,v*3.6:9:1,v:6:1);
    showsound;
    ht:=ht+0.5;
   end;
 end
else showsound;

end.