c===========================================================
c     democomplex: Demonstration of basic complex number 
c     manipulations in Fortran 77.
c
c     Sample usage on lnx1:
c
c     lnx1% pwd
c     /home/phys410/util
c
c     lnx1% make democomplex
c     pgf77 -g -L/usr/local/PGI/lib \
c        democomplex.o -o democomplex
c     Linking:
c
c     lnx1% democomplex
c      i =  (0.0000000000000000E+000,1.000000000000000)
c      a =  (3.000000000000000,3.000000000000000)
c      b =  (3.000000000000000,0.0000000000000000E+000)
c      i**2 =  (-1.000000000000000,0.0000000000000000E+000)
c      a / b =  (1.000000000000000,1.000000000000000)
c      6.0d0 * a * b - i =  (54.00000000000000,53.00000000000000)
c      c =  (2.000000000000000,3.000000000000000)
c      dble(i) =   0.0000000000000000E+000
c      imag(i) =    1.000000000000000
c      abs(i) =    1.000000000000000
c      abs(a) =    4.242640687119285
c      conjg(a) =  (3.000000000000000,-3.000000000000000
c
c===========================================================
      program    democomplex

      implicit   none

c-----------------------------------------------------------
c     complex*16 is the equivalent of real*8, i.e. real*8
c     values are used for both the real and imaginary parts
c     of the complex numbers.  Complex arrays can be defined 
c     as usual.
c-----------------------------------------------------------
      complex*16 a, b, i, c, ary1(10), ary2(10,10), 
     &           ary3(10,10,10)
      real*8     realc, imagc

c-----------------------------------------------------------
c     complex*16 constants are specified as 
c 
c     ( <real-part> , <imaginary-part> )
c-----------------------------------------------------------
      i = (0.0d0,1.0d0)
      a = (3.0d0,3.0d0)
      b = (3.0d0,0.0d0)

      write(*,*) 'i = ', i
      write(*,*) 'a = ', a
      write(*,*) 'b = ', b

c-----------------------------------------------------------
c     Basic arithmetic on complex numbers is fully 
c     supported.
c-----------------------------------------------------------
      write(*,*) 'i**2 = ', i**2
      write(*,*) 'a / b = ', a / b
      write(*,*) '6.0d0 * a * b - i = ', 6.0d0 * a * b - i

c-----------------------------------------------------------
c     Use the 'dcmplx' intrinsic to convert real and imaginary
c     parts into a single complex*16 value.
c-----------------------------------------------------------
      realc = 2.0d0
      imagc = 3.0d0
      c     = dcmplx(realc,imagc)
      write(*,*) 'c = ', c

c-----------------------------------------------------------
c     Intrinsic functions for manipulations of complex 
c     numbers:
c
c     dble(z)  returns   the real part of z
c     imag(z)  returns   the imaginary part of z
c     abs(z)   returns   the modulus of z, i.e. sqrt(z z*)
c     conjg(z) returns   the complex conjugate of z
c-----------------------------------------------------------
      write(*,*) 'dble(i) = ', dble(i)
      write(*,*) 'imag(i) = ', imag(i)
      write(*,*) 'abs(i) = ', abs(i)
      write(*,*) 'abs(a) = ', abs(a)
      write(*,*) 'conjg(a) = ', conjg(a)

      stop

      end
