DataMuseum.dk

boolean procedure err int cox (x, y, u, v);
_______ _________
  value                        x, y       ;
  _____
  real                         x, y, u, v ;
  ____
    comment  err int cox  computes the value of the
    _______
    ERRor INTegral with COmpleX argument.
    The parameters are:
    x:           real      part of input z ,
    y:           imaginary part of input z ,
    u:           real      part of output w ,
    v:           imaginary part of output w ,
    err int cox: is true when 0 _ arg(z) _ phi/4 ,
                    ____        <        <
                 otherwise it is false;
                                 _____
if x < 0 ∨ y < 0 ∨ y > 1.00001⨯x
__
  then err int cox := false
  ____                _____
  else
  ____
    begin
    _____
      err int cox := true;
                     ____
      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;
            ____
            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;
              ___
            u := y⨯a;
            v := b
          end Hermite quadrature
          ___
        else
        ____
          begin comment Legendre approximation;
          _____ _______
            real p1,p2,p3,n1,n2,t1,t2,a,b,T,M;
            ____
            procedure PK(pa,pb,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10);
            _________
              value            a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 ;
              _____
              real       pa,pb,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 :=  12096.51250
                        + a ⨯ p1  + M ⨯ p2   ;
              pb :=       b ⨯ p1
            end PK;
            ___
            a := (x + y)⨯(x - y);
            b := 2⨯x⨯y;
            T := 0.4⨯a;
            M := - 0.04⨯(a∧2 + b∧2);
                          |     |
            PK(t1   ,t2   ,  -8488.78070,  14448.00988,
               -4495.93759,   3287.20821,   -519.3045 ,
                 210.21   ,    -14.3    ,      3.3    ,
                   0      ,      0                     );
            PK(n1   ,n2   ,  31832.92763,  39914.35198,
               31537.26576,  17481.0636 ,   7151.3442 ,
                2207.205  ,    514.8    ,     89.1    ,
                  11      ,      1                     );
            p3 := 1.12837 91671 /(n1∧2 + n2∧2);
                                    |      |
            p2 := p3⨯(n1⨯t1 + n2⨯t2);
            p1 := p3⨯(n1⨯t2 - n2⨯t1);
            T  := exp(-a);
            u  :=   T⨯cos(b) - x⨯p1 - y⨯p2;
            v  := - T⨯sin(b) + x⨯p2 - y⨯p1
          end Legendre approximation
          ___
    end 0 _ arg(z) _ phi/4
    ___   <        <
finis   err int cox;
[end] [stop]