DataMuseum.dk

boolean procedure Sommerfeld cox (pr, pi, gr, gi);
_______ _________
  value                           pr, pi         ;
  _____
  real                            pr, pi, gr, gi ;
  ____
    comment This procedure computes the value of
    _______
    the Sommerfeld attenuation function: G(p).
    The parameters are:
    pr: real      part of input p,
    pi: imaginary part of input p,
    gr: real      part of output G(p),
    gi: imaginary part of output G(p),
    Sommerfeld cox: is true when 0 _ arg(p) _ phi/2,
                       ____        <        <
                    otherwise it is false;
                                    _____
if pr < 0 ∨ pi < 0
__
  then Sommerfeld cox := false
  ____                   _____
  else
  ____
    begin
    _____
      real x,y,M;
      ____
      Sommerfeld cox := true;
                        ____
      M := pr∧2 + pi∧2;
             |      |
      x := sqrt((sqrt(M) + pr)/2);
      y := if x = 0 then 0 else pi/2/x;
           __       ____   ____
      if y > 1.7 - 0.2 ⨯ x ∨ y > 3.9 - x
      __
        then
        ____
          begin comment Hermite quadrature;
          _____ _______
            real p1,p2,p3,p4,p5,p6,n1,n2,n3,n4,n5,n6,a,b,T;
            ____
            M := y∧2;
                  |
            a := b := 0;
            for T := - x , x do
            ___              __
              begin
              _____
                p1 := 0.31424 03763 + T     ;
                p2 := 0.94778 83912 + T     ;
                p3 := 1.59768 26352 + T     ;
                p4 := 2.27950 70805 + T     ;
                p5 := 3.02063 70251 + T     ;
                p6 := 3.88972 48979 + T     ;
                n1 := 0.18147 96822         /(p1∧2 + M);
                                                |
                n2 := 0.08291 72776 3       /(p2∧2 + M);
                                                |
                n3 := 0.01642 73320 3       /(p3∧2 + M);
                                                |
                n4 := 0.00124 31244 32      /(p4∧2 + M);
                                                |
                n5 := 0.00002 72908 9347    /(p5∧2 + M);
                                                |
                n6 := 0.00000 00846 24328 41/(p6∧2 + M);
                                                |
                a  :=   a + n1 + n2 + n3 + n4 + n5 + n6;
                b  := - b + p1⨯n1 + p2⨯n2 + p3⨯n3 +
                            p4⨯n4 + p5⨯n5 + p6⨯n6
              end T;
              ___
            gr := 1 - 1.77245 38509⨯(x⨯b + M⨯a);
            gi :=     1.77245 38509⨯(x⨯a -   b)⨯y
          end Hermite quadrature
          ___
        else
        ____
          begin comment Legendre approximation;
          _____ _______
            real p1,p2,p3,n1,n2,t1,t2,T;
            ____
            procedure PK(pa,pb,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10);
            _________
              value            a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 ;
              _____
              real       pa,pb,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 ;
              ____
            begin
            _____
              p3 :=  a9 + T ⨯ a10            ;
              p2 :=  a8 + T ⨯ p3  + M ⨯ a10  ;
              p1 :=  a7 + T ⨯ p2  + M ⨯ p3   ;
              p3 :=  a6 + T ⨯ p1  + M ⨯ p2   ;
              p2 :=  a5 + T ⨯ p3  + M ⨯ p1   ;
              p1 :=  a4 + T ⨯ p2  + M ⨯ p3   ;
              p3 :=  a3 + T ⨯ p1  + M ⨯ p2   ;
              p2 :=  a2 + T ⨯ p3  + M ⨯ p1   ;
              p1 := (a1 + T ⨯ p2  + M ⨯ p3)/5;
              pa :=  a0 + pr⨯ p1  + M ⨯ p2   ;
              pb :=       pi⨯ p1
            end PK;
            ___
            T := 0.4⨯pr;
            M := - 0.04⨯M;
            PK(t1   ,t2   ,
                   0      ,  12096.51250,  -8488.78070,
               14448.00988,  -4495.93759,   3287.20821,
                -519.3045 ,    210.21   ,    -14.3    ,
                   3.3    ,      0                     );
            PK(n1   ,n2   ,
               12096.51250,  31832.92763,  39914.35198,
               31537.26576,  17481.0636 ,   7151.3442 ,
                2207.205  ,    514.8    ,     89.1    ,
                  11      ,      1                     );
            p3 := 10/(n1∧2 + n2∧2);
                        |      |
            p2 := cos(pi);
            p1 := sin(pi);
            T  := 1.77245 38509⨯exp(-pr);
            gr := 1 + T⨯(x⨯p1 - y⨯p2) - p3⨯(n1⨯t1 + n2⨯t2);
            gi :=     T⨯(x⨯p2 + y⨯p1) - p3⨯(n1⨯t2 - n2⨯t1)
          end Legendre approximation
          ___
    end 0 _ arg(p) _ phi/2
    ___   <        <
finis   Sommerfeld cox;
[end] [stop]