home *** CD-ROM | disk | FTP | other *** search
/ Whiteline: Alpha / Whiteline Alpha.iso / progtool / modula2 / complex / cmplx_16.for next >
Encoding:
Text File  |  1994-09-22  |  9.9 KB  |  384 lines
  1.  
  2. ! Bibliothek von Unterprogrammen für doppelt-genaue komplexe Rechnungen
  3. ! Die Beschreibung finden Sie in CMPLX_16.TXT . Zerrupfen Sie die Samm-
  4. ! lung nur mit Bedacht, da die Programme z.T. aufeinander zugreifen.
  5. !       (C) 1989  Thomas Groß, Lynarstraße 13, D 1 Berlin 65
  6.  
  7.  
  8.       SUBROUTINE CDACOS(A,B,C,D)
  9.  
  10.       REAL*8 A,B,C,D
  11.       CALL CDARCUS(A,B,C,D)
  12.       C=ACOS(C)
  13.       D=-D
  14.       END
  15.  
  16. !---------------------------------------------------------------------!
  17.  
  18.       SUBROUTINE CDACOSH(A,B,C,D)
  19.  
  20.       REAL*8 A,B,C,D
  21.       CALL CDACOS(A,B,D,C)                          ! Abramowitz 4.6.15
  22.       C=-C
  23.       END
  24.  
  25. !---------------------------------------------------------------------!
  26.  
  27.       SUBROUTINE CDACOT(A,B,C,D)
  28.  
  29.       REAL*8 A,B,C,D
  30.       CALL CDATAN(-A,-B,C,D)                         ! Abramowitz 4.4.5
  31.       C=IVORZEI(A)*ASIN(1.D0)+C
  32.       END
  33.  
  34. !---------------------------------------------------------------------!
  35.  
  36.       SUBROUTINE CDACOTH(A,B,C,D)
  37.  
  38.       REAL*8 A,B,C,D
  39.       CALL CDACOT(-B,A,D,C)                         ! Abramowitz 4.6.19
  40.       C=-C
  41.       END
  42.  
  43. !---------------------------------------------------------------------!
  44.  
  45.       SUBROUTINE CDASIN(A,B,C,D)
  46.  
  47.       REAL*8 A,B,C,D
  48.       CALL CDARCUS(A,B,C,D)
  49.       C=ASIN(C)
  50.       END
  51.  
  52. !---------------------------------------------------------------------!
  53.  
  54.       SUBROUTINE CDASINH(A,B,C,D)
  55.  
  56.       REAL*8 A,B,C,D
  57.       CALL CDASIN(-B,A,D,C)                         ! Abramowitz 4.6.14
  58.       D=-D
  59.       END
  60.  
  61. !---------------------------------------------------------------------!
  62.  
  63.       SUBROUTINE CDATAN(A,B,C,D)
  64.  
  65.       REAL*8 A,B,C,D
  66.       IF(B) 2,1,2
  67.     1 C=ATAN(A)
  68.       D=0.
  69.       RETURN
  70.     2 D=1-A*A-B*B
  71.       IF(D) 7,3,7
  72.     3 IF(A) 6,4,6
  73.     4 C=0.
  74.     5 D=IVORZEI(B)*1.D308
  75.       RETURN
  76.     6 C=IVORZEI(A)*ASIN(1.D0)/2
  77.       IF(ABS(A).GE.1.5D-154 .OR. ABS(B)-1.NE.0.) GOTO 8
  78.       PRINT *
  79.       PRINT *,'        Aus dem Unterprogramm CDATAN: Underflow von LOG(A
  80.      &^2) !'
  81.       PAUSE'- Die Rückgabewerte sind  Re = sgn(A) π/4  und  Im = sgn(B) 
  82.      &1.D308'
  83.       GOTO 5      
  84.     7 C=ATAN2(2*A,D)/2                              ! Abramowitz 4.4.39
  85.     8 D=LOG((A*A+(B+1)**2)/(A*A+(B-1)**2))/4
  86.       END
  87.  
  88. !---------------------------------------------------------------------!
  89.  
  90.       SUBROUTINE CDATANH(A,B,C,D)
  91.  
  92.       REAL*8 A,B,C,D
  93.       CALL CDATAN(-B,A,D,C)                         ! Abramowitz 4.6.16
  94.       D=-D
  95.       END
  96.  
  97. !---------------------------------------------------------------------!
  98.  
  99.       SUBROUTINE CDCOS(A,B,C,D)
  100.  
  101.       REAL*8 A,B,C,D,CD0D0,CD7D2
  102.       C=CD0D0(COS(A))*COSH(CD7D2(B))
  103.       D=SIN(A)
  104.       IF(ABS(A).GE.1.) D=CD0D0(D)
  105.       D=-D*SINH(CD7D2(B))
  106.       END
  107.  
  108. !---------------------------------------------------------------------!
  109.       
  110.       SUBROUTINE CDCOSH(A,B,C,D)
  111.  
  112.       REAL*8 A,B,C,D
  113.       CALL CDCOS(-B,A,C,D)                           ! Abramowitz 4.5.8
  114.       END
  115.  
  116. !---------------------------------------------------------------------!
  117.       
  118.       SUBROUTINE CDCOT(A,B,C,D)
  119.  
  120.       REAL*8 A,B,C,D,V,W,X,Y,Z,CD0D0
  121.       IF(B.NE.0.) GOTO 4
  122.       D=0.
  123.       Y=COS(A)
  124.       Z=SIN(A)
  125.       IF(ABS(A).GE.1.) Z=CD0D0(Z)
  126.       IF(Z.NE.0.) GOTO 3
  127.       IF(Y) 1,2,2
  128.     1 C=-1.D308
  129.       RETURN
  130.     2 C=1.D308
  131.       RETURN
  132.     3 C=CD0D0(Y)/Z
  133.       RETURN
  134.     4 V=B
  135.       IF(ABS(V).GT.355.) V=IVORZEI(V)*3.55D2   ! kein Overflow in CDDIV
  136.       CALL CDCOS(A,V,W,X)
  137.       CALL CDSIN(A,V,Y,Z)
  138.       CALL CDDIV(W,X,Y,Z,C,D)
  139.       END
  140.  
  141. !---------------------------------------------------------------------!
  142.       
  143.       SUBROUTINE CDCOTH(A,B,C,D)
  144.  
  145.       REAL*8 A,B,C,D
  146.       CALL CDCOT(-B,A,D,C)                          ! Abramowitz 4.5.12
  147.       C=-C
  148.       END
  149.       
  150. !---------------------------------------------------------------------!
  151.       
  152.       SUBROUTINE CDDIV(A1,B1,A2,B2,C,D)
  153.  
  154.       REAL*8 A1,A2,B1,B2,C,D
  155.       D=A2*A2+B2*B2
  156.       IF(A2.NE.0. .OR. B2.NE.0.) GOTO 2
  157.       PRINT *
  158.       PRINT *,'        Aus dem Unterprogramm CDDIV: Der Nenner ist 0 .'
  159.     1 PAUSE'- Die Rückgabewerte sind  Re = 0.  und  Im = 0.'
  160.       C=0.
  161.       RETURN
  162.     2 IF(D) 4,3,4
  163.     3 PRINT *
  164.       PRINT *,'        Aus dem Unterprogramm CDDIV: Underflow des Nenner
  165.      &betrags!'
  166.       GOTO 1
  167.     4 C=(A1*A2+B1*B2)/D
  168.       D=(A2*B1-A1*B2)/D       
  169.       END
  170.  
  171. !---------------------------------------------------------------------!
  172.  
  173.       SUBROUTINE CDEXP(A,B,C,D)
  174.  
  175.       REAL*8 A,B,C,D,CD0D0,CD7D2
  176.       C=EXP(CD7D2(A))
  177.       D=SIN(B)
  178.       IF(ABS(B).GE.1.) D=CD0D0(D)
  179.       D=C*D
  180.       C=C*CD0D0(COS(B))
  181.       END
  182.  
  183. !---------------------------------------------------------------------!
  184.       
  185.       SUBROUTINE CDLOG(A,B,C,D)
  186.  
  187.       REAL*8 A,B,C,D
  188.       IF(A.NE.0.) GOTO 5
  189.       IF(B) 1,3,2
  190.     1 C=LOG(-B)
  191.       D=-ASIN(1.D0)
  192.       RETURN
  193.     2 C=LOG(B)
  194.       D=ASIN(1.D0)
  195.       RETURN
  196.     3 PRINT *
  197.       PRINT *,'        Aus dem Unterprogramm CDLOG:  ln 0  ist nicht def
  198.      &iniert!'
  199.     4 PAUSE'- Die Rückgabewerte sind  Re = - 1.D308  und  Im = 0.'
  200.       C=-1.D308
  201.       D=0.
  202.       RETURN
  203.     5 C=SQRT(A*A+B*B)
  204.       IF(C) 7,6,7
  205.     6 PRINT *
  206.       PRINT *,'        Aus dem Unterprogramm CDLOG: Underflow von SQRT(A
  207.      &^2 + B^2) !'
  208.       GOTO 4
  209.     7 C=LOG(C)                                       ! Abramowitz 4.1.2
  210.       D=ATAN2(B,A)
  211.       END
  212.  
  213. !---------------------------------------------------------------------!
  214.       
  215.       SUBROUTINE CDMUL(A1,B1,A2,B2,C,D)
  216.  
  217.       REAL*8 A1,A2,B1,B2,C,D
  218.       C=A1*A2-B1*B2
  219.       D=A1*B2+A2*B1
  220.       END
  221.  
  222. !---------------------------------------------------------------------!
  223.       
  224.       SUBROUTINE CDPOT(A1,B1,A2,B2,C,D)
  225.  
  226.       REAL*8 A1,A2,B1,B2,C,D
  227.       IF(A1.NE.0..OR.B1.NE.0.) GOTO 2
  228.       IF(A2.NE.0..OR.B2.NE.0.) GOTO 1
  229.       PRINT *
  230.       PRINT *,'        Aus dem Unterprogramm CDPOT:  0^0  ist nicht defi
  231.      &niert!'
  232.       PAUSE'- Die Rückgabewerte sind  Re = 0.  und  Im = 0.'
  233.     1 C=0.
  234.       D=0.
  235.       RETURN
  236.     2 CALL CDLOG(A1,B1,C,D)                          ! Abramowitz 4.2.7
  237.       CALL CDEXP(A2*C-B2*D,A2*D+B2*C,C,D)
  238.       END
  239.  
  240. !---------------------------------------------------------------------!
  241.       
  242.       SUBROUTINE CDSIN(A,B,C,D)
  243.  
  244.       REAL*8 A,B,C,D,CD0D0,CD7D2
  245.       C=SIN(A)
  246.       IF(ABS(A).GE.1.) C=CD0D0(C)
  247.       C=C*COSH(CD7D2(B))
  248.       D=CD0D0(COS(A))*SINH(CD7D2(B))
  249.       END
  250.  
  251. !---------------------------------------------------------------------!
  252.       
  253.       SUBROUTINE CDSINH(A,B,C,D)
  254.  
  255.       REAL*8 A,B,C,D
  256.       CALL CDSIN(-B,A,D,C)                           ! Abramowitz 4.5.7
  257.       D=-D
  258.       END
  259.       
  260. !---------------------------------------------------------------------!
  261.       
  262.       SUBROUTINE CDSQRT(A,B,C,D)
  263.  
  264.       REAL*8 A,B,C,D,Z,CD0D0
  265.       IF(A.NE.0.) GOTO 3
  266.       C=SQRT(ABS(B)/2)
  267.       IF(B) 1,2,2
  268.     1 D=-C
  269.       RETURN
  270.     2 D=C
  271.       RETURN
  272.     3 D=SQRT(SQRT(A*A+B*B))                         ! Abramowitz 3.7.26
  273.       IF(D) 5,4,5
  274.     4 PRINT *
  275.       PRINT *,'        Aus dem Unterprogramm CDSQRT: Underflow von SQRT(
  276.      &A^2 + B^2) !'
  277.       PAUSE'- Die Rückgabewerte sind  Re = 0.  und  Im = 0.'
  278.     5 Z=ATAN2(B,A)/2
  279.       C=D*CD0D0(COS(Z))
  280.       D=D*SIN(Z)
  281.       END
  282.  
  283. !---------------------------------------------------------------------!
  284.       
  285.       SUBROUTINE CDTAN(A,B,C,D)
  286.  
  287.       REAL*8 A,B,C,D,V,W,X,Y,Z,CD0D0
  288.       IF(B.NE.0.) GOTO 1
  289.       C=TAN(A)
  290.       IF(ABS(A).GE.1.) C=CD0D0(C)
  291.       IF(ABS(C).GT.1.E16) C=IVORZEI(C)*1.D308     ! Vergrößert auf  ± ∞
  292.       D=0.
  293.       RETURN
  294.     1 V=B
  295.       IF(ABS(V).GT.355.) V=IVORZEI(V)*3.55D2   ! kein Overflow in CDDIV
  296.       CALL CDSIN(A,V,W,X)
  297.       CALL CDCOS(A,V,Y,Z)
  298.       CALL CDDIV(W,X,Y,Z,C,D)
  299.       END
  300.  
  301. !---------------------------------------------------------------------!
  302.       
  303.       SUBROUTINE CDTANH(A,B,C,D)
  304.  
  305.       REAL*8 A,B,C,D
  306.       CALL CDTAN(-B,A,D,C)                           ! Abramowitz 4.5.9
  307.       D=-D
  308.       END
  309.  
  310. !=====================================================================!
  311.  
  312. ! Hilfsprogramme
  313.  
  314.  
  315.       REAL*8 FUNCTION CD0D0(X)
  316. ! Verkleinert  < ± 2.D-16  auf  0.
  317. ! wegen Rundungsfehlern von cos(π/2) , sin(π) , tan(π)
  318.       REAL*8 X
  319.       CD0D0=X
  320.       IF(ABS(X).LT.2.E-16) CD0D0=0.
  321.       END
  322.  
  323.  
  324.       REAL*8 FUNCTION CD7D2(X)
  325. ! Verkleinert  > ± 704  auf  ± 704  wegen Overflow
  326. ! von exp( ) , cosh( ) , sinh( )
  327.       REAL*8 X
  328.       CD7D2=X
  329.       IF(ABS(X).GT.704.) CD7D2=IVORZEI(X)*7.04D2
  330.       END
  331.  
  332.  
  333.       SUBROUTINE CDARCUS(A,B,Z,D)
  334.       REAL*8 A,B,C,D,Y,Z,CDKETTE
  335.       C=SQRT((1+A)**2+B*B)                       ! Abramowitz 4.4.37,38
  336.       D=SQRT((1-A)**2+B*B)
  337.       Y=(C+D)/2
  338.       Z=A
  339.       IF(B) 2,1,2
  340.     1 IF(ABS(A).LT.1.D0) GOTO 6
  341.       Y=ABS(A)
  342.       Z=IVORZEI(A)
  343.       GOTO 6
  344.     2 IF(A) 4,3,4
  345.     3 D=IVORZEI(B)*LOG(Y+ABS(B))
  346.       RETURN
  347.     4 Z=(C-D)/2      ! numerisch sehr gefährliche Differenz der Wurzeln
  348.       V=ABS(B/A)
  349.       IF(ABS(A).LE.2. .OR. V.GT.0.3 .AND. V.LT.3.) GOTO 5
  350.       IF(V.LE.0.3) THEN
  351.        C=1/(A+1)                   ! Reihenentwicklung für 0,3|A| > |B|
  352.        D=1/(A-1)
  353.        Z=B*B/4
  354.        Z=IVORZEI(A)*(1-ABS(Z*(C*CDKETTE(C*Z*C)-D*CDKETTE(D*Z*D))))
  355.       ELSE
  356.        C=(A+1)**2                  ! Reihenentwicklung für 3|A| < |B|
  357.        D=(A-1)**2
  358.        Z=4*B*B
  359.        Z=(C*CDKETTE(C/Z)-D*CDKETTE(D/Z))/4/ABS(B)
  360.       ENDIF
  361.     5 IF(ABS(Z).GT.1.D0) Z=IVORZEI(A)            ! gegen Rundungsfehler
  362.     6 D=IVORZEI(B)*LOG(Y+SQRT(Y*Y-1))
  363.       END
  364.  
  365.  
  366.       REAL*8 FUNCTION CDKETTE(X)
  367. ! Reihenentwicklung der Wurzel in CDARCUS , Abramowitz 3.6.11
  368.       REAL*8 X
  369.       CDKETTE=1+X*(-1+X*(2+X*(-5+X*(14+X*(-42+X*(132+X*(-429+
  370.      &          X*(1430+X*(-4862+X*(16796-X*48155))))))))))       
  371.       END         ! eigene Restgliedkorrektur ^
  372.  
  373.  
  374.       FUNCTION IVORZEI(X)
  375. ! spezielle Vorzeichenfunktion
  376.       REAL*8 X
  377.       IF(X) 1,2,2
  378.     1 IVORZEI=-1
  379.       RETURN
  380.     2 IVORZEI=1
  381.       END
  382.  
  383. !                               * * *
  384.