home *** CD-ROM | disk | FTP | other *** search
/ PC World 2000 August / PCWorld_2000-08_cd.bin / Software / TemaCD / xbasic / xbpro.exe / xb / xcm.x < prev    next >
Text File  |  2000-01-12  |  25KB  |  996 lines
  1. '
  2. '
  3. ' ####################  Max Reason
  4. ' #####  PROLOG  #####  copyright 1988-2000
  5. ' ####################  XBasic complex number function library
  6. '
  7. ' subject to GPLL license - see gpll.txt
  8. '
  9. ' maxreason@maxreason.com
  10. '
  11. ' for Windows XBasic
  12. ' for Linux XBasic
  13. '
  14. PROGRAM    "xcm"
  15. VERSION    "0.0007"
  16. '
  17. IMPORT    "xst"
  18. IMPORT    "xma"
  19. '
  20. EXPORT
  21. '
  22. ' *******************************
  23. ' *****  declare functions  *****
  24. ' *******************************
  25. '
  26. DECLARE FUNCTION Xcm                   ()
  27. DECLARE FUNCTION XcmVersion$           ()
  28. '
  29. ' DCOMPLEX functions  (.R and .I components are DOUBLE precision)
  30. '
  31. DECLARE FUNCTION DOUBLE    DCABS       (DCOMPLEX z)
  32. DECLARE FUNCTION DCOMPLEX  DCACOS      (DCOMPLEX z)
  33. DECLARE FUNCTION DOUBLE    DCARG       (DCOMPLEX z)
  34. DECLARE FUNCTION DCOMPLEX  DCASIN      (DCOMPLEX z)
  35. DECLARE FUNCTION DCOMPLEX  DCATAN      (DCOMPLEX z)
  36. DECLARE FUNCTION DCOMPLEX  DCCONJ      (DCOMPLEX z)
  37. DECLARE FUNCTION DCOMPLEX  DCCOS       (DCOMPLEX z)
  38. DECLARE FUNCTION DCOMPLEX  DCCOSH      (DCOMPLEX z)
  39. DECLARE FUNCTION DCOMPLEX  DCEXP       (DCOMPLEX z)
  40. DECLARE FUNCTION DCOMPLEX  DCLOG       (DCOMPLEX z)
  41. DECLARE FUNCTION DCOMPLEX  DCLOG10     (DCOMPLEX z)
  42. DECLARE FUNCTION DOUBLE    DCNORM      (DCOMPLEX z)
  43. DECLARE FUNCTION DCOMPLEX  DCPOLAR     (DOUBLE mag, DOUBLE angle)
  44. DECLARE FUNCTION DCOMPLEX  DCPOWERCC   (DCOMPLEX z, DCOMPLEX n)
  45. DECLARE FUNCTION DCOMPLEX  DCPOWERCR   (DCOMPLEX z, DOUBLE   n)
  46. DECLARE FUNCTION DCOMPLEX  DCPOWERRC   (DOUBLE   z, DCOMPLEX n)
  47. DECLARE FUNCTION DCOMPLEX  DCRMUL      (DCOMPLEX x, DOUBLE y)
  48. DECLARE FUNCTION DCOMPLEX  DCSIN       (DCOMPLEX z)
  49. DECLARE FUNCTION DCOMPLEX  DCSINH      (DCOMPLEX z)
  50. DECLARE FUNCTION DCOMPLEX  DCSQRT      (DCOMPLEX z)
  51. DECLARE FUNCTION DCOMPLEX  DCTAN       (DCOMPLEX z)
  52. DECLARE FUNCTION DCOMPLEX  DCTANH      (DCOMPLEX z)
  53. '
  54. ' SCOMPLEX functions  (.R and .I components are SINGLE precision)
  55. '
  56. DECLARE FUNCTION SINGLE    SCABS       (SCOMPLEX z)
  57. DECLARE FUNCTION SCOMPLEX  SCACOS      (SCOMPLEX z)
  58. DECLARE FUNCTION SINGLE    SCARG       (SCOMPLEX z)
  59. DECLARE FUNCTION SCOMPLEX  SCASIN      (SCOMPLEX z)
  60. DECLARE FUNCTION SCOMPLEX  SCATAN      (SCOMPLEX z)
  61. DECLARE FUNCTION SCOMPLEX  SCCONJ      (SCOMPLEX z)
  62. DECLARE FUNCTION SCOMPLEX  SCCOS       (SCOMPLEX z)
  63. DECLARE FUNCTION SCOMPLEX  SCCOSH      (SCOMPLEX z)
  64. DECLARE FUNCTION SCOMPLEX  SCEXP       (SCOMPLEX z)
  65. DECLARE FUNCTION SCOMPLEX  SCLOG       (SCOMPLEX z)
  66. DECLARE FUNCTION SCOMPLEX  SCLOG10     (SCOMPLEX z)
  67. DECLARE FUNCTION SINGLE    SCNORM      (SCOMPLEX z)
  68. DECLARE FUNCTION SCOMPLEX  SCPOLAR     (SINGLE mag, SINGLE angle)
  69. DECLARE FUNCTION SCOMPLEX  SCPOWERCC   (SCOMPLEX z, SCOMPLEX n)
  70. DECLARE FUNCTION SCOMPLEX  SCPOWERCR   (SCOMPLEX z, SINGLE   n)
  71. DECLARE FUNCTION SCOMPLEX  SCPOWERRC   (SINGLE   z, SCOMPLEX n)
  72. DECLARE FUNCTION SCOMPLEX  SCRMUL      (SCOMPLEX x, SINGLE y)
  73. DECLARE FUNCTION SCOMPLEX  SCSIN       (SCOMPLEX z)
  74. DECLARE FUNCTION SCOMPLEX  SCSINH      (SCOMPLEX z)
  75. DECLARE FUNCTION SCOMPLEX  SCSQRT      (SCOMPLEX z)
  76. DECLARE FUNCTION SCOMPLEX  SCTAN       (SCOMPLEX z)
  77. DECLARE FUNCTION SCOMPLEX  SCTANH      (SCOMPLEX z)
  78. '
  79. END EXPORT
  80. '
  81. INTERNAL FUNCTION DOUBLE  XdcGetAlpha    (DCOMPLEX z)      ' called by DCASIN & DCACOS
  82. INTERNAL FUNCTION DOUBLE  XdcGetBeta     (DCOMPLEX z)      ' ditto
  83. INTERNAL FUNCTION SINGLE  XscGetAlpha    (SCOMPLEX z)      ' called by SCASIN & SCACOS
  84. INTERNAL FUNCTION SINGLE  XscGetBeta     (SCOMPLEX z)      ' ditto
  85. INTERNAL FUNCTION DOUBLE  Atan2          (DOUBLE y, DOUBLE x)
  86. '
  87. '
  88. ' ####################
  89. ' #####  Xcm ()  #####
  90. ' ####################
  91. '
  92. FUNCTION  Xcm () DOUBLE
  93.     XLONG i, j, ofile
  94.     DCOMPLEX z, zneg, zexp
  95.     DCOMPLEX polar
  96.     DCOMPLEX conj, sq1, sq2, sq3
  97.     DCOMPLEX expx, expxi, emeix, epeix, log1, log10x
  98.     DCOMPLEX sinr, cosr, tanr, asinx, acosx, atanx
  99.     DCOMPLEX sinhx, coshx, tanhx
  100.     DCOMPLEX powCC, powCR, powRC
  101. '
  102.     RETURN
  103. '
  104.     a$ = "Max Reason"
  105.     a$ = "copyright 1988-2000"
  106.     a$ = "XBasic complex number function library"
  107.     a$ = "maxreason@maxreason.com"
  108.     a$ = ""
  109. '
  110. ' Comment out the following line to route output to file ofile$
  111. '
  112.     ofile = $$STDOUT
  113.     IF (ofile != $$STDOUT) THEN
  114.         ofile$ = "complex.out"
  115.          ofile = OPEN (ofile$, $$WRNEW)
  116.       IF (ofile <= 0) THEN PRINT [$$STDERR], "cannot open ofile"
  117.     END IF
  118. '
  119.     DO
  120.       a$ = INLINE$ ("Input base real value      ===>> ")
  121.         IFZ a$ THEN EXIT FUNCTION
  122.         b$ = INLINE$ ("Input base imaginary value ===>> ")
  123.         z.R = DOUBLE (a$)
  124.         z.I = DOUBLE (b$)
  125. '
  126. ' Uncomment the following 4 lines when testing DCPOWERDC, DCPOWERCR or DCPOWERRC
  127. '
  128.       a$ = INLINE$ ("Input exp real value       ===>> ")
  129.         b$ = INLINE$ ("Input exp imaginary value  ===>> ")
  130.         zexp.R = DOUBLE (a$)
  131.         zexp.I = DOUBLE (b$)
  132.  
  133.         zneg.R = -z.R
  134.         zneg.I = -z.I
  135.  
  136.         norm   = DCNORM(z)
  137.         aabs         = DCABS(z)
  138.         arg    = DCARG(z)
  139.         polar  = DCPOLAR(aabs, arg)
  140.         conj     = DCCONJ(z)
  141.         sq1         = DCSQRT(z)
  142.         log10x = DCLOG10(z)
  143.         log1   = DCLOG(z)
  144.         expx   = DCEXP(z)
  145.         expxi  = DCEXP(zneg)
  146.         emeix  = expx - expxi
  147.         epeix  = expx + expxi
  148.         sinr     = DCSIN(z)
  149.         cosr   = DCCOS(z)
  150.         tanr   = DCTAN(z)
  151.         asinx     = DCASIN(sinr)
  152.         acosx  = DCACOS(cosr)
  153.         atanx  = DCATAN(tanr)
  154.         asinx     = DCASIN(z)
  155.         acosx  = DCACOS(z)
  156.         atanx  = DCATAN(z)
  157.         sinhx  = DCSINH(z)
  158.         coshx  = DCCOSH(z)
  159.         tanhx  = DCTANH(z)
  160.  
  161.         powCC  = DCPOWERCC(z, zexp)
  162.         powCR  = DCPOWERCR(z, zexp.R)
  163.         powRC  = DCPOWERRC(z.R, zexp)
  164.  
  165.         PRINT [ofile], "abs      = "; aabs,  TAB(36); HEX$(DHIGH(aabs), 8);;  HEX$(DLOW(aabs), 8)
  166.         PRINT [ofile], "arg      = "; arg,   TAB(36); HEX$(DHIGH(arg), 8);;   HEX$(DLOW(arg), 8)
  167.         PRINT [ofile], "norm     = "; norm,  TAB(36); HEX$(DHIGH(norm), 8);;  HEX$(DLOW(norm), 8)
  168.         PRINT [ofile], "polar.R  = "; polar.R,  TAB(36); HEX$(DHIGH(polar.R), 8);;  HEX$(DLOW(polar.R), 8)
  169.         PRINT [ofile], "polar.I  = "; polar.I,  TAB(36); HEX$(DHIGH(polar.I), 8);;  HEX$(DLOW(polar.I), 8)
  170.         PRINT [ofile], "conj.R   = "; conj.R,  TAB(36); HEX$(DHIGH(conj.R), 8);;  HEX$(DLOW(conj.R), 8)
  171.         PRINT [ofile], "conj.I   = "; conj.I,  TAB(36); HEX$(DHIGH(conj.I), 8);;  HEX$(DLOW(conj.I), 8)
  172.         PRINT [ofile], "sq1.R    = "; sq1.R,  TAB(36); HEX$(DHIGH(sq1.R), 8);;  HEX$(DLOW(sq1.R), 8)
  173.         PRINT [ofile], "sq1.I    = "; sq1.I,  TAB(36); HEX$(DHIGH(sq1.I), 8);;  HEX$(DLOW(sq1.I), 8)
  174.         PRINT [ofile], "log1.R   = "; log1.R,  TAB(36); HEX$(DHIGH(log1.R), 8);;  HEX$(DLOW(log1.R), 8)
  175.         PRINT [ofile], "log1.I   = "; log1.I,  TAB(36); HEX$(DHIGH(log1.I), 8);;  HEX$(DLOW(log1.I), 8)
  176.         PRINT [ofile], "log10x.R = "; log10x.R,  TAB(36); HEX$(DHIGH(log10x.R), 8);;  HEX$(DLOW(log10x.R), 8)
  177.         PRINT [ofile], "log10x.I = "; log10x.I,  TAB(36); HEX$(DHIGH(log10x.I), 8);;  HEX$(DLOW(log10x.I), 8)
  178.         PRINT [ofile], "expx.R   = "; expx.R,  TAB(36); HEX$(DHIGH(expx.R), 8);;  HEX$(DLOW(expx.R), 8)
  179.         PRINT [ofile], "expx.I   = "; expx.I,  TAB(36); HEX$(DHIGH(expx.I), 8);;  HEX$(DLOW(expx.I), 8)
  180.         PRINT [ofile], "expxi.R  = "; expxi.R,  TAB(36); HEX$(DHIGH(expxi.R), 8);;  HEX$(DLOW(expxi.R), 8)
  181.         PRINT [ofile], "expxi.I  = "; expxi.I,  TAB(36); HEX$(DHIGH(expxi.I), 8);;  HEX$(DLOW(expxi.I), 8)
  182.         PRINT [ofile], "emeix.R  = "; emeix.R,  TAB(36); HEX$(DHIGH(emeix.R), 8);;  HEX$(DLOW(emeix.R), 8)
  183.         PRINT [ofile], "emeix.I  = "; emeix.I,  TAB(36); HEX$(DHIGH(emeix.I), 8);;  HEX$(DLOW(emeix.I), 8)
  184.         PRINT [ofile], "epeix.R  = "; epeix.R,  TAB(36); HEX$(DHIGH(epeix.R), 8);;  HEX$(DLOW(epeix.R), 8)
  185.         PRINT [ofile], "epeix.I  = "; epeix.I,  TAB(36); HEX$(DHIGH(epeix.I), 8);;  HEX$(DLOW(epeix.I), 8)
  186.         PRINT [ofile], "sinr.R   = "; sinr.R,  TAB(36); HEX$(DHIGH(sinr.R), 8);;  HEX$(DLOW(sinr.R), 8)
  187.         PRINT [ofile], "sinr.I   = "; sinr.I,  TAB(36); HEX$(DHIGH(sinr.I), 8);;  HEX$(DLOW(sinr.I), 8)
  188.         PRINT [ofile], "cosr.R   = "; cosr.R,  TAB(36); HEX$(DHIGH(cosr.R), 8);;  HEX$(DLOW(cosr.R), 8)
  189.         PRINT [ofile], "cosr.I   = "; cosr.I,  TAB(36); HEX$(DHIGH(cosr.I), 8);;  HEX$(DLOW(cosr.I), 8)
  190.         PRINT [ofile], "tanr.R   = "; tanr.R,  TAB(36); HEX$(DHIGH(tanr.R), 8);;  HEX$(DLOW(tanr.R), 8)
  191.         PRINT [ofile], "tanr.I   = "; tanr.I,  TAB(36); HEX$(DHIGH(tanr.I), 8);;  HEX$(DLOW(tanr.I), 8)
  192.         PRINT [ofile], "asin.R   = "; asinx.R, TAB(36); HEX$(DHIGH(asinx.R), 8);; HEX$(DLOW(asinx.R), 8)
  193.         PRINT [ofile], "asin.I   = "; asinx.I, TAB(36); HEX$(DHIGH(asinx.I), 8);; HEX$(DLOW(asinx.I), 8)
  194.         PRINT [ofile], "acos.R   = "; acosx.R, TAB(36); HEX$(DHIGH(acosx.R), 8);; HEX$(DLOW(acosx.R), 8)
  195.         PRINT [ofile], "acos.I   = "; acosx.I, TAB(36); HEX$(DHIGH(acosx.I), 8);; HEX$(DLOW(acosx.I), 8)
  196.         PRINT [ofile], "atan.R   = "; atanx.R, TAB(36); HEX$(DHIGH(atanx.R), 8);; HEX$(DLOW(atanx.R), 8)
  197.         PRINT [ofile], "atan.I   = "; atanx.I, TAB(36); HEX$(DHIGH(atanx.I), 8);; HEX$(DLOW(atanx.I), 8)
  198.         PRINT [ofile], "sinhx.R  = "; sinhx.R,  TAB(36); HEX$(DHIGH(sinhx.R), 8);;  HEX$(DLOW(sinhx.R), 8)
  199.         PRINT [ofile], "sinhx.I  = "; sinhx.I,  TAB(36); HEX$(DHIGH(sinhx.I), 8);;  HEX$(DLOW(sinhx.I), 8)
  200.         PRINT [ofile], "coshx.R  = "; coshx.R,  TAB(36); HEX$(DHIGH(coshx.R), 8);;  HEX$(DLOW(coshx.R), 8)
  201.         PRINT [ofile], "coshx.I  = "; coshx.I,  TAB(36); HEX$(DHIGH(coshx.I), 8);;  HEX$(DLOW(coshx.I), 8)
  202.         PRINT [ofile], "tanhx.R  = "; tanhx.R,  TAB(36); HEX$(DHIGH(tanhx.R), 8);;  HEX$(DLOW(tanhx.R), 8)
  203.         PRINT [ofile], "tanhx.I  = "; tanhx.I,  TAB(36); HEX$(DHIGH(tanhx.I), 8);;  HEX$(DLOW(tanhx.I), 8)
  204.  
  205.         PRINT [ofile], "powCC.R  = "; powCC.R, TAB(36); HEX$(DHIGH(powCC.R), 8);; HEX$(DLOW(powCC.R), 8)
  206.         PRINT [ofile], "powCC.I  = "; powCC.I, TAB(36); HEX$(DHIGH(powCC.I), 8);; HEX$(DLOW(powCC.I), 8)
  207.         PRINT [ofile], "powCR.R  = "; powCR.R,  TAB(36); HEX$(DHIGH(powCR.R), 8);;  HEX$(DLOW(powCR.R), 8)
  208.         PRINT [ofile], "powCR.I  = "; powCR.I,  TAB(36); HEX$(DHIGH(powCR.I), 8);;  HEX$(DLOW(powCR.I), 8)
  209.         PRINT [ofile], "powRC.R  = "; powRC.R,  TAB(36); HEX$(DHIGH(powRC.R), 8);;  HEX$(DLOW(powRC.R), 8)
  210.         PRINT [ofile], "powRC.I  = "; powRC.I,  TAB(36); HEX$(DHIGH(powRC.I), 8);;  HEX$(DLOW(powRC.I), 8)
  211.     LOOP
  212.   IF (ofile != $$STDOUT) THEN CLOSE (ofile)
  213. END FUNCTION
  214. '
  215. '
  216. '  ############################
  217. '  #####  XcmVersion$ ()  #####
  218. '  ############################
  219. '
  220. FUNCTION  XcmVersion$ ()
  221.     version$ = VERSION$ (0)
  222.     RETURN (version$)
  223. END FUNCTION
  224. '
  225. '
  226. ' ######################
  227. ' #####  DCABS ()  ##### this version uses the Numerical Recipes method
  228. ' ######################
  229. '
  230. FUNCTION DOUBLE DCABS (DCOMPLEX z) DOUBLE
  231. '
  232.     x = ABS(z.R)
  233.     y = ABS(z.I)
  234. '
  235.     SELECT CASE TRUE
  236.         CASE (x = 0#):    RETURN y
  237.         CASE (y = 0#):    RETURN x
  238.         CASE (x > y):        tmp = y / x : RETURN (x * SQRT(1#+(tmp*tmp)))
  239.         CASE ELSE:            tmp = x / y : RETURN (y * SQRT(1#+(tmp*tmp)))
  240.     END SELECT
  241. END FUNCTION
  242. '
  243. '
  244. ' #######################
  245. ' #####  DCACOS ()  ##### from Handbook of Mathematical Functions
  246. ' #######################
  247. '
  248. FUNCTION DCOMPLEX DCACOS (DCOMPLEX z) DOUBLE
  249.     DCOMPLEX ans
  250. '
  251.     alpha = XdcGetAlpha(z)
  252.     beta  = XdcGetBeta (z)
  253.     ans.R = ACOS(beta)
  254.     ans.I = -alpha
  255.     RETURN (ans)
  256. END FUNCTION
  257. '
  258. '
  259. ' ######################
  260. ' #####  DCARG ()  ##### returns radians
  261. ' ######################
  262. '
  263. FUNCTION DOUBLE DCARG (DCOMPLEX z) DOUBLE
  264.     IF (z.I = 0) & (z.R = 0) THEN RETURN (0#)
  265.     RETURN (Atan2(z.I, z.R))
  266. END FUNCTION
  267. '
  268. '
  269. ' #######################
  270. ' #####  DCASIN ()  ##### from Handbook of Math Functions
  271. ' #######################
  272. '
  273. FUNCTION DCOMPLEX DCASIN (DCOMPLEX z) DOUBLE
  274.     DCOMPLEX ans
  275. '
  276.     alpha = XdcGetAlpha(z)
  277.     beta  = XdcGetBeta (z)
  278.  
  279.     ans.R = ASIN(beta)
  280.     ans.I = alpha
  281.     RETURN (ans)
  282. END FUNCTION
  283. '
  284. '
  285. ' #######################
  286. ' #####  DCATAN ()  #####
  287. ' #######################
  288. '
  289. FUNCTION DCOMPLEX DCATAN (DCOMPLEX z) DOUBLE
  290.     DCOMPLEX    ans
  291.     XLONG     ##ERROR
  292. '
  293. '    z*z may not = -1"
  294. '
  295.     IF ((z.R = 0#) & (z.I = 1#)) THEN ##ERROR = $$ErrorNatureInvalidArgument: RETURN (0#)
  296. '
  297.     xsq = z.R * z.R
  298.     ysq = z.I * z.I
  299.     yincsq = z.I + 1#
  300.     yincsq = yincsq * yincsq
  301.     ydecsq = z.I - 1#
  302.     ydecsq = ydecsq * ydecsq
  303.     ans.R = ATAN((z.R * 2#) / (1# - xsq - ysq)) * 0.5#
  304.     ans.I = LOG((xsq + yincsq) / (xsq + ydecsq)) * 0.25#
  305.     RETURN (ans)
  306. END FUNCTION
  307. '
  308. '
  309. ' #######################
  310. ' #####  DCCONJ ()  #####
  311. ' #######################
  312. '
  313. FUNCTION DCOMPLEX DCCONJ (DCOMPLEX z) DOUBLE
  314.     DCOMPLEX ans
  315. '
  316.     ans.R =  z.R
  317.     ans.I = -z.I
  318.     RETURN (ans)
  319. END FUNCTION
  320. '
  321. '
  322. ' ######################
  323. ' #####  DCCOS ()  #####
  324. ' ######################
  325. '
  326. FUNCTION DCOMPLEX DCCOS (DCOMPLEX z) DOUBLE
  327.     DCOMPLEX ans
  328. '
  329.     ans.R =  COS(z.R) * COSH(z.I)
  330.     ans.I = -SIN(z.R) * SINH(z.I)
  331.     RETURN (ans)
  332. END FUNCTION
  333. '
  334. '
  335. ' #######################
  336. ' #####  DCCOSH ()  #####
  337. ' #######################
  338. '
  339. FUNCTION DCOMPLEX DCCOSH (DCOMPLEX z) DOUBLE
  340.     DCOMPLEX ans
  341. '
  342.     ans.R = COSH(z.R) * COS(z.I)
  343.     ans.I = SINH(z.R) * SIN(z.I)
  344.     RETURN (ans)
  345. END FUNCTION
  346. '
  347. '
  348. ' ######################
  349. ' #####  DCEXP ()  #####
  350. ' ######################
  351. '
  352. FUNCTION DCOMPLEX DCEXP (DCOMPLEX z) DOUBLE
  353.     DCOMPLEX ans
  354. '
  355.     ans.R = EXP(z.R) * COS(z.I)
  356.     ans.I = EXP(z.R) * SIN(z.I)
  357.     RETURN (ans)
  358. END FUNCTION
  359. '
  360. '
  361. ' ######################
  362. ' #####  DCLOG ()  #####
  363. ' ######################
  364. '
  365. FUNCTION DCOMPLEX DCLOG (DCOMPLEX z) DOUBLE
  366.     DCOMPLEX ans
  367. '
  368.     ans.R = LOG(DCNORM(z)) * 0.5#
  369.     ans.I = Atan2(z.I, z.R)
  370.     RETURN (ans)
  371. END FUNCTION
  372. '
  373. '
  374. ' ########################
  375. ' #####  DCLOG10 ()  #####
  376. ' ########################
  377. '
  378. FUNCTION DCOMPLEX DCLOG10 (DCOMPLEX z) DOUBLE
  379.     DCOMPLEX ans
  380. '
  381.     ans = DCLOG(z)
  382.     ans.R = ans.R * $$LOG10E
  383.     ans.I = ans.I * $$LOG10E
  384.     RETURN (ans)
  385. END FUNCTION
  386. '
  387. '
  388. ' #######################
  389. ' #####  DCNORM ()  #####
  390. ' #######################
  391. '
  392. FUNCTION DOUBLE DCNORM (DCOMPLEX z) DOUBLE
  393. '
  394.     RETURN ((z.R * z.R) + (z.I * z.I))
  395. END FUNCTION
  396. '
  397. '
  398. ' ########################
  399. ' #####  DCPOLAR ()  #####
  400. ' ########################
  401. '
  402. FUNCTION DCOMPLEX DCPOLAR (DOUBLE mag, DOUBLE angle) DOUBLE
  403.     DCOMPLEX ans
  404. '
  405.     ans.R = mag * COS(angle)
  406.     ans.I = mag * SIN(angle)
  407.     RETURN (ans)
  408. END FUNCTION
  409. '
  410. '
  411. ' ##########################
  412. ' #####  DCPOWERCC ()  #####
  413. ' ##########################
  414. '
  415. FUNCTION DCOMPLEX DCPOWERCC (DCOMPLEX z, DCOMPLEX n) DOUBLE
  416.     DCOMPLEX ans
  417. '
  418.     ans = DCEXP( n * DCLOG(z) )
  419.     RETURN (ans)
  420. END FUNCTION
  421. '
  422. '
  423. ' ##########################
  424. ' #####  DCPOWERCR ()  #####
  425. ' ##########################
  426. '
  427. FUNCTION DCOMPLEX DCPOWERCR (DCOMPLEX z, DOUBLE n) DOUBLE
  428.     DCOMPLEX ans
  429. '
  430.     ans = DCEXP( DCRMUL(DCLOG(z), n) )
  431.     RETURN (ans)
  432. END FUNCTION
  433. '
  434. '
  435. ' ##########################
  436. ' #####  DCPOWERRC ()  #####
  437. ' ##########################
  438. '
  439. FUNCTION DCOMPLEX DCPOWERRC (DOUBLE z, DCOMPLEX n) DOUBLE
  440.     DCOMPLEX ans
  441. '
  442.     ans = DCEXP( DCRMUL(n, LOG(z)) )
  443.     RETURN (ans)
  444. END FUNCTION
  445. '
  446. '
  447. ' #######################
  448. ' #####  DCRMUL ()  #####
  449. ' #######################
  450. '
  451. FUNCTION DCOMPLEX DCRMUL (DCOMPLEX x, DOUBLE y) DOUBLE
  452.     DCOMPLEX ans
  453. '
  454.     ans.R = x.R * y
  455.     ans.I = x.I * y
  456.     RETURN (ans)
  457. END FUNCTION
  458. '
  459. '
  460. ' ######################
  461. ' #####  DCSIN ()  #####
  462. ' ######################
  463. '
  464. FUNCTION DCOMPLEX DCSIN (DCOMPLEX z) DOUBLE
  465.     DCOMPLEX ans
  466. '
  467.     ans.R = SIN(z.R) * COSH(z.I)
  468.     ans.I = COS(z.R) * SINH(z.I)
  469.     RETURN (ans)
  470. END FUNCTION
  471. '
  472. '
  473. ' #######################
  474. ' #####  DCSINH ()  #####
  475. ' #######################
  476. '
  477. FUNCTION DCOMPLEX DCSINH (DCOMPLEX z) DOUBLE
  478.     DCOMPLEX ans
  479. '
  480.     ans.R = SINH(z.R) * COS(z.I)
  481.     ans.I = COSH(z.R) * SIN(z.I)
  482.     RETURN (ans)
  483. END FUNCTION
  484. '
  485. '
  486. ' #######################
  487. ' #####  DCSQRT ()  #####  Math Handbook
  488. ' #######################
  489. '
  490. FUNCTION DCOMPLEX DCSQRT (DCOMPLEX z) DOUBLE
  491.     DCOMPLEX ans
  492. '
  493.     zabs    = DCABS(z)
  494.     ans.R    = SQRT( (zabs + z.R) * 0.5#)
  495.     ans.I    = SQRT( (zabs - z.R) * 0.5#) * SIGN(z.I)
  496.     RETURN (ans)
  497. '
  498. '    *****************************************
  499. '    *****  DCSQRT2:  numeric recipes:  *****
  500. '    *****************************************
  501. '
  502. '    IF (z.R = 0#) & (z.I = 0#) THEN ans.R = 0# : ans.I = 0# : RETURN (ans)
  503. '    x = ABS(z.R)
  504. '    y = ABS(z.I)
  505. '    IF (x >= y) THEN
  506. '        r = y/x
  507. '        w = SQRT(x) * SQRT(0.5# * (1# + SQRT(1# + (r*r))))
  508. '    ELSE
  509. '        r = x/y
  510. '        w = SQRT(y) * SQRT(0.5# * (r  + SQRT(1# + (r*r))))
  511. '    END IF
  512. '    IF z.R >= 0 THEN
  513. '        ans.R = w
  514. '        ans.I = z.I / (2# * w)
  515. '    ELSE
  516. '        IF z.I >= 0 THEN ans.I = w ELSE ans.I = -w
  517. '        ans.R = z.I / (2# * ans.I)
  518. '    END IF
  519. '    RETURN (ans)
  520. '
  521. '    ****************************************
  522. '    *****  DCSQRT3:  Turbo C++ method  *****
  523. '    ****************************************
  524. '
  525. '    sqrtzabs = SQRT(DCABS(z))
  526. '    hargz = DCARG(z) / 2#
  527. '    ans.R = sqrtzabs * COS(hargz)
  528. '    ans.I = sqrtzabs * SIN(hargz)
  529. '    RETURN (ans)
  530. '
  531. END FUNCTION
  532. '
  533. '
  534. ' ######################
  535. ' #####  DCTAN ()  #####
  536. ' ######################
  537. '
  538. FUNCTION DCOMPLEX DCTAN (DCOMPLEX z) DOUBLE
  539.     DCOMPLEX ans
  540. '
  541.     denom = 1# / (COS(2# * z.R) + COSH(2# * z.I))
  542.     ans.R = SIN (2# * z.R) * denom
  543.     ans.I = SINH (2# * z.I) * denom
  544.     RETURN (ans)
  545. END FUNCTION
  546. '
  547. '
  548. ' #######################
  549. ' #####  DCTANH ()  #####
  550. ' #######################
  551. '
  552. FUNCTION DCOMPLEX DCTANH (DCOMPLEX z) DOUBLE
  553.     DCOMPLEX ans
  554. '
  555.     denom = 1# / (COSH(2# * z.R) + COS(2# * z.I))
  556.     ans.R = SINH(2# * z.R) * denom
  557.     ans.I = SIN (2# * z.I) * denom
  558.     RETURN (ans)
  559. END FUNCTION
  560. '
  561. '
  562. ' ######################
  563. ' #####  SCABS ()  #####  This version uses the Numerical Recipes method
  564. ' ######################
  565. '
  566. FUNCTION SINGLE SCABS (SCOMPLEX z) SINGLE
  567. '
  568.     x = ABS(z.R)
  569.     y = ABS(z.I)
  570. '
  571.     SELECT CASE TRUE
  572.         CASE (x = 0#):    RETURN y
  573.         CASE (y = 0#):    RETURN x
  574.         CASE (x > y):        tmp = y / x : RETURN (x * SQRT(1#+(tmp*tmp)))
  575.         CASE ELSE:            tmp = x / y : RETURN (y * SQRT(1#+(tmp*tmp)))
  576.     END SELECT
  577. END FUNCTION
  578. '
  579. '
  580. ' #######################
  581. ' #####  SCACOS ()  ##### from Handbook of Mathematical Functions
  582. ' #######################
  583. '
  584. FUNCTION SCOMPLEX SCACOS (SCOMPLEX z) SINGLE
  585.     SCOMPLEX ans
  586. '
  587.     alpha = XscGetAlpha(z)
  588.     beta  = XscGetBeta (z)
  589.     ans.R = ACOS(beta)
  590.     ans.I = -alpha
  591.     RETURN (ans)
  592. END FUNCTION
  593. '
  594. '
  595. ' ######################
  596. ' #####  SCARG ()  ##### returns radians
  597. ' ######################
  598. '
  599. FUNCTION SINGLE SCARG (SCOMPLEX z) SINGLE
  600. '
  601.     IF (z.I = 0) & (z.R = 0) THEN RETURN (0#)
  602.     RETURN (Atan2(z.I, z.R))
  603. END FUNCTION
  604. '
  605. '
  606. ' #######################
  607. ' #####  SCASIN ()  ##### from Handbook of Math Functions
  608. ' #######################
  609. '
  610. FUNCTION SCOMPLEX SCASIN (SCOMPLEX z) SINGLE
  611.     SCOMPLEX ans
  612. '
  613.     alpha = XscGetAlpha(z)
  614.     beta  = XscGetBeta (z)
  615.     ans.R = ASIN(beta)
  616.     ans.I = alpha
  617.     RETURN (ans)
  618. END FUNCTION
  619. '
  620. '
  621. ' #######################
  622. ' #####  SCATAN ()  #####
  623. ' #######################
  624. '
  625. FUNCTION SCOMPLEX SCATAN (SCOMPLEX z) SINGLE
  626.     SCOMPLEX    ans
  627.     XLONG     ##ERROR
  628. '
  629. '    z*z may not = -1"
  630. '
  631.     IF ((z.R = 0#) & (z.I = 1#)) THEN
  632.         ##ERROR = $$ErrorNatureInvalidArgument
  633.         RETURN (0#)
  634.     END IF
  635. '
  636.     xsq = z.R * z.R
  637.     ysq = z.I * z.I
  638. '
  639.     yincsq = z.I + 1#
  640.     yincsq = yincsq * yincsq
  641. '
  642.     ydecsq = z.I - 1#
  643.     ydecsq = ydecsq * ydecsq
  644. '
  645.     ans.R = ATAN((z.R * 2#) / (1# - xsq - ysq)) * 0.5#
  646.     ans.I = LOG((xsq + yincsq) / (xsq + ydecsq)) * 0.25#
  647.     RETURN (ans)
  648. END FUNCTION
  649. '
  650. '
  651. ' #######################
  652. ' #####  SCCONJ ()  #####
  653. ' #######################
  654. '
  655. FUNCTION SCOMPLEX SCCONJ (SCOMPLEX z) SINGLE
  656.     SCOMPLEX ans
  657. '
  658.     ans.R =     z.R
  659.     ans.I = - z.I
  660.     RETURN (ans)
  661. END FUNCTION
  662. '
  663. '
  664. ' ######################
  665. ' #####  SCCOS ()  #####
  666. ' ######################
  667. '
  668. FUNCTION SCOMPLEX SCCOS (SCOMPLEX z) SINGLE
  669.     SCOMPLEX ans
  670. '
  671.     ans.R =  COS(z.R) * COSH(z.I)
  672.     ans.I = -SIN(z.R) * SINH(z.I)
  673.     RETURN (ans)
  674. END FUNCTION
  675. '
  676. '
  677. ' #######################
  678. ' #####  SCCOSH ()  #####
  679. ' #######################
  680. '
  681. FUNCTION SCOMPLEX SCCOSH (SCOMPLEX z) SINGLE
  682.     SCOMPLEX ans
  683. '
  684.     ans.R = COSH(z.R) * COS(z.I)
  685.     ans.I = SINH(z.R) * SIN(z.I)
  686.     RETURN (ans)
  687. END FUNCTION
  688. '
  689. '
  690. ' ######################
  691. ' #####  SCEXP ()  #####
  692. ' ######################
  693. '
  694. FUNCTION SCOMPLEX SCEXP (SCOMPLEX z) SINGLE
  695.     SCOMPLEX ans
  696. '
  697.     ans.R = EXP(z.R) * COS(z.I)
  698.     ans.I = EXP(z.R) * SIN(z.I)
  699.     RETURN (ans)
  700. END FUNCTION
  701. '
  702. '
  703. ' ######################
  704. ' #####  SCLOG ()  #####
  705. ' ######################
  706. '
  707. FUNCTION SCOMPLEX SCLOG (SCOMPLEX z) SINGLE
  708.     SCOMPLEX ans
  709. '
  710.     ans.R = LOG(SCNORM(z)) * 0.5#
  711.     ans.I = Atan2(z.I, z.R)
  712.     RETURN (ans)
  713. END FUNCTION
  714. '
  715. '
  716. ' ########################
  717. ' #####  SCLOG10 ()  #####
  718. ' ########################
  719. '
  720. FUNCTION SCOMPLEX SCLOG10 (SCOMPLEX z) SINGLE
  721.     SCOMPLEX ans
  722. '
  723.     ans = SCLOG(z)
  724.     ans.R = ans.R * $$LOG10E
  725.     ans.I = ans.I * $$LOG10E
  726.     RETURN (ans)
  727. END FUNCTION
  728. '
  729. '
  730. ' #######################
  731. ' #####  SCNORM ()  #####
  732. ' #######################
  733. '
  734. FUNCTION SINGLE SCNORM (SCOMPLEX z) SINGLE
  735.     RETURN ((z.R * z.R) + (z.I * z.I))
  736. END FUNCTION
  737. '
  738. '
  739. ' ########################
  740. ' #####  SCPOLAR ()  #####
  741. ' ########################
  742. '
  743. FUNCTION SCOMPLEX SCPOLAR (SINGLE mag, SINGLE angle) SINGLE
  744.     SCOMPLEX ans
  745. '
  746.     ans.R = mag * COS(angle)
  747.     ans.I = mag * SIN(angle)
  748.     RETURN (ans)
  749. END FUNCTION
  750. '
  751. '
  752. ' ##########################
  753. ' #####  SCPOWERCC ()  #####
  754. ' ##########################
  755. '
  756. FUNCTION SCOMPLEX SCPOWERCC (SCOMPLEX z, SCOMPLEX n) SINGLE
  757.     SCOMPLEX ans
  758. '
  759.     ans = SCEXP( n * SCLOG(z) )
  760.     RETURN (ans)
  761. END FUNCTION
  762. '
  763. '
  764. ' ##########################
  765. ' #####  SCPOWERCR ()  #####
  766. ' ##########################
  767. '
  768. FUNCTION SCOMPLEX SCPOWERCR (SCOMPLEX z, SINGLE n) SINGLE
  769.     SCOMPLEX ans
  770. '
  771.     ans = SCEXP( SCRMUL(SCLOG(z), n) )
  772.     RETURN (ans)
  773. END FUNCTION
  774. '
  775. '
  776. ' ##########################
  777. ' #####  SCPOWERRC ()  #####
  778. ' ##########################
  779. '
  780. FUNCTION SCOMPLEX SCPOWERRC (SINGLE z, SCOMPLEX n) SINGLE
  781.     SCOMPLEX ans
  782. '
  783.     ans = SCEXP( SCRMUL(n, LOG(z)) )
  784.     RETURN (ans)
  785. END FUNCTION
  786. '
  787. '
  788. ' #######################
  789. ' #####  SCRMUL ()  #####
  790. ' #######################
  791. '
  792. FUNCTION SCOMPLEX SCRMUL (SCOMPLEX x, SINGLE y) SINGLE
  793.     SCOMPLEX ans
  794. '
  795.     ans.R = x.R * y
  796.     ans.I = x.I * y
  797.     RETURN (ans)
  798. END FUNCTION
  799. '
  800. '
  801. ' ######################
  802. ' #####  SCSIN ()  #####
  803. ' ######################
  804. '
  805. FUNCTION SCOMPLEX SCSIN (SCOMPLEX z) SINGLE
  806.     SCOMPLEX ans
  807. '
  808.     ans.R = SIN(z.R) * COSH(z.I)
  809.     ans.I = COS(z.R) * SINH(z.I)
  810.     RETURN (ans)
  811. END FUNCTION
  812. '
  813. '
  814. ' #######################
  815. ' #####  SCSINH ()  #####
  816. ' #######################
  817. '
  818. FUNCTION SCOMPLEX SCSINH (SCOMPLEX z) SINGLE
  819.     SCOMPLEX ans
  820. '
  821.     ans.R = SINH(z.R) * COS(z.I)
  822.     ans.I = COSH(z.R) * SIN(z.I)
  823.     RETURN (ans)
  824. END FUNCTION
  825. '
  826. '
  827. ' #######################
  828. ' #####  SCSQRT ()  #####  Math Handbook
  829. ' #######################
  830. '
  831. FUNCTION SCOMPLEX SCSQRT (SCOMPLEX z) SINGLE
  832.     SCOMPLEX ans
  833. '
  834.     zabs    = SCABS(z)
  835.     ans.R    = SQRT( (zabs + z.R) * 0.5#)
  836.     ans.I    = SQRT( (zabs - z.R) * 0.5#) * SIGN(z.I)
  837.     RETURN (ans)
  838. '
  839. '    *****************************************
  840. '    *****  XSqrt2:  numeric recipes:  *****
  841. '    *****************************************
  842. '
  843. '    IF (z.R = 0#) & (z.I = 0#) THEN ans.R = 0# : ans.I = 0# : RETURN (ans)
  844. '    x = ABS(z.R)
  845. '    y = ABS(z.I)
  846. '    IF (x >= y) THEN
  847. '        r = y/x
  848. '        w = SQRT(x) * SQRT(0.5# * (1# + SQRT(1# + (r*r))))
  849. '    ELSE
  850. '        r = x/y
  851. '        w = SQRT(y) * SQRT(0.5# * (r  + SQRT(1# + (r*r))))
  852. '    END IF
  853. '    IF z.R >= 0 THEN
  854. '        ans.R = w
  855. '        ans.I = z.I / (2# * w)
  856. '    ELSE
  857. '        IF z.I >= 0 THEN ans.I = w ELSE ans.I = -w
  858. '        ans.R = z.I / (2# * ans.I)
  859. '    END IF
  860. '    RETURN (ans)
  861. '
  862. '    ****************************************
  863. '    *****  SCSQRT3:  Turbo C++ method  *****
  864. '    ****************************************
  865. '
  866. '    sqrtzabs = SQRT(SCABS(z))
  867. '    hargz = SCARG(z) / 2#
  868. '    ans.R = sqrtzabs * COS (hargz)
  869. '    ans.I = sqrtzabs * SIN (hargz)
  870. '    RETURN (ans)
  871. '
  872. END FUNCTION
  873. '
  874. '
  875. ' ######################
  876. ' #####  SCTAN ()  #####
  877. ' ######################
  878. '
  879. FUNCTION SCOMPLEX SCTAN (SCOMPLEX z) SINGLE
  880.     SCOMPLEX ans
  881. '
  882.     denom = 1# / (COS(2# * z.R) + COSH(2# * z.I))
  883.     ans.R = SIN (2# * z.R) * denom
  884.     ans.I = SINH(2# * z.I) * denom
  885.     RETURN (ans)
  886. END FUNCTION
  887. '
  888. '
  889. ' #######################
  890. ' #####  SCTANH ()  #####
  891. ' #######################
  892. '
  893. FUNCTION SCOMPLEX SCTANH (SCOMPLEX z) SINGLE
  894.     SCOMPLEX ans
  895. '
  896.     denom = 1# / (COSH(2# * z.R) + COS(2# * z.I))
  897.     ans.R = SINH(2# * z.R) * denom
  898.     ans.I = SIN (2# * z.I) * denom
  899.     RETURN (ans)
  900. END FUNCTION
  901. '
  902. '
  903. ' ############################  From Handbook of Mathematical Functions
  904. ' #####  XdcGetAlpha ()  #####  Called by XAsin and DCACOS
  905. ' ############################
  906. '
  907. '    returns the complete imaginary term, to be specific
  908. '
  909. FUNCTION DOUBLE XdcGetAlpha (DCOMPLEX z) DOUBLE
  910.     ysq = z.I * z.I
  911.     xincsq = z.R + 1#
  912.     xincsq = xincsq * xincsq
  913.     xdecsq = z.R - 1#
  914.     xdecsq = xdecsq * xdecsq
  915.     alpha = (SQRT( xincsq + ysq ) + SQRT( xdecsq + ysq ) ) * 0.5#
  916.     RETURN ( LOG(alpha + SQRT((alpha*alpha) - 1#)) )
  917. END FUNCTION
  918. '
  919. '
  920. ' ###########################
  921. ' #####  XdcGetBeta ()  #####
  922. ' ###########################
  923. '
  924. FUNCTION DOUBLE XdcGetBeta (DCOMPLEX z) DOUBLE
  925.     ysq = z.I * z.I
  926.     xincsq = z.R + 1#
  927.     xincsq = xincsq * xincsq
  928.     xdecsq = z.R - 1#
  929.     xdecsq = xdecsq * xdecsq
  930.     beta = (SQRT( xincsq + ysq ) - SQRT( xdecsq + ysq ) ) * 0.5#
  931.     RETURN (beta)
  932. END FUNCTION
  933. '
  934. '
  935. ' ############################  From Handbook of Mathematical Functions
  936. ' #####  XscGetAlpha ()  #####  Called by SCASIN and XAcos
  937. ' ############################
  938. '
  939. '    returns the complete imaginary term, to be specific
  940. '
  941. FUNCTION SINGLE XscGetAlpha (SCOMPLEX z) SINGLE
  942.     ysq = z.I * z.I
  943.     xincsq = z.R + 1#
  944.     xincsq = xincsq * xincsq
  945.     xdecsq = z.R - 1#
  946.     xdecsq = xdecsq * xdecsq
  947.     alpha = (SQRT( xincsq + ysq ) + SQRT( xdecsq + ysq ) ) * 0.5#
  948.     RETURN ( LOG(alpha + SQRT((alpha*alpha) - 1#)) )
  949. END FUNCTION
  950. '
  951. '
  952. ' ###########################
  953. ' #####  XscGetBeta ()  #####
  954. ' ###########################
  955. '
  956. FUNCTION SINGLE XscGetBeta (SCOMPLEX z) SINGLE
  957.     ysq = z.I * z.I
  958.     xincsq = z.R + 1#
  959.     xincsq = xincsq * xincsq
  960.     xdecsq = z.R - 1#
  961.     xdecsq = xdecsq * xdecsq
  962.     beta = (SQRT( xincsq + ysq ) - SQRT( xdecsq + ysq ) ) * 0.5#
  963.     RETURN (beta)
  964. END FUNCTION
  965. '
  966. '
  967. ' ######################
  968. ' #####  Atan2 ()  #####
  969. ' ######################
  970. '
  971. FUNCTION DOUBLE Atan2 (DOUBLE y, DOUBLE x) DOUBLE
  972.     XLONG     ##ERROR
  973. '
  974.     IF ((x = 0) & (y = 0)) THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN(0#)
  975.     IF ((y = 0) & (x < 0)) THEN RETURN ($$PI)
  976. '
  977.     xsign = SIGN(x)
  978.     ysign = SIGN(y)
  979.  
  980.     x = ABS(x)
  981.     y = ABS(y)
  982.  
  983.     IF (x >= y) THEN
  984.         s = (y/x)
  985.     ELSE
  986.         s = (x/y)
  987.         inv = 1
  988.     END IF
  989.  
  990.     res = ATAN(s)
  991.     IF inv THEN res = $$PIDIV2 - res            ' atan(y/x) = pi/2 - atan(x/y)
  992.     IF (xsign = -1#) THEN res = $$PI - res
  993.     RETURN (ysign * res)
  994. END FUNCTION
  995. END PROGRAM
  996.