Complex bakeoff - example 3

[Jim Thomas]

Complex bakeoff - example 3                              10/12/95     3:27 PM

Bakeoff example 1 used complex infinity.  Example 2 used multiple 
complex infinities.  This example uses the sign of zero.

In "Augmenting a Programming Language with Complex Arithmetic" (NCEG/91-039), 
W. Kahan and I presented two examples illustrating how the sign of zero 
can be used to distinguish the two sides of a branch cut to simplify the 
conceptualization and coding for problems involving branch cuts.  Below
is a simple, straightforward C program for the first of the two examples, 
but the program depends on consistent treatment of the sign of zero.  The
paper includes plots showing adverse affects of treating -0 like +0.  

This example is intended to illustrate a basic use of the sign of zero.  More 
complicated examples better dramatize the relative difficulty in developing 
a solution that does not utilize the sign of zero.  Still, for comparison 
purposes, how might one code this examples for the Cray proposal which 
does not support the sign of zero?  Does the code become more complicated 
or less efficient?  What is the effort to devise the workaround?

In "Handling of Complex Numbers in the C^H Programming Language" 
(X3J11.1/93-020), section 7.7, Harry H. Cheng addresses these very examples, 
with reference to his experimental extended C implementation.  He rightly 
concludes that distinguishing the signs of zero is not necessary to program 
the two examples.  However, his schemes are more complicated, and obviously 
required time to devise and prove.

Example 3

As a conformal map,  w = f(z) = (z - 1/z) / 2  maps the complex z-plane 
twice, once for  |z| >= 1  and once for   |z| <= 1 , onto the complex w-plane,

mapping the unit circle  |z| = 1  to a slit along the imaginary axis from  
w = -i  to w = +i.  The inverse map

        g(w) = w - i*sqrt(i*w - 1)*sqrt(i*w + 1)

maps the whole w-plane slitted along  w = -i  to w = +i  onto the outside of
the 
z-plane's unit circle  |z| >= 1 .  Vertical lines in the w-plane map to the 
streamlines of a vertical "eluding" flow around the unit circle on the
z-plane.
We wish to calculate points to exhibit these streamlines.

To smooth the plot near stagnation points  g(+/-i)  where  g'(+/-i)  is
infinite
we use

        h(s) = s*(s*s*(3*s*s - 10) + 15) / 8

which satisfies  h(-s) = -h(s) ,  h'(s) >= 0 ,  h(1) = 1 ,  and h'(1) = 0 .
For any fixed real  r ,  w = r + i*h(s)  traces out a vertical straight line 
segment as we let  s  run through  -1.5 <= s <= 1.5  , in small steps
ds = 3/32.  z = g(w)  traces out a streamline past the circle.
        
To plot several streamlines, we run  r  through  0 <= r <= 0.6 , in steps of 
dr = 0.05 .  For each such  r  we plot two streamlines,  z = g(r + i*h(s))
to the right of the circle and  z = g(-r + i*h(s))  to the left.

The main program calculates the points  z  on the streamlines and prints 
the real and imaginary parts into two files.  MATLAB can be used to load 
the file of imaginary parts and the file of real parts into array variables, 
reshape them into  (2*N+1)x(2*M+2)  arrays  X  and  Y , respectively, and 
then plot column of  Y   against column of  X .

//-------------------------------------------------------
//      Eluding flow past a disk
//      using a C++ prototype of "Complex C Extensions"
//-------------------------------------------------------

#include <stdio.h>
#include <complex.h>                    // Complex C Extensions

// Constants:

const int       M       = 12;           // Number of streamlines = 2*M+1 (one
splits) 
const int       N       = 16;           // Number of points per steamline =
2*N+1
const double    dr      = 0.05;         // Distance between streamlines
const double    ds      = 3.0 / 32.0;   // Step-size

// Functions:

double_complex g(double_complex w) { return w - I * sqrt(I * w - 1) * sqrt(I *
w + 1); }
double h(double s) { return s * (s * s * (3 * s * s - 10) + 15) / 8 ; }
void print_parts(double_complex z);

// To generate points along streamlines ...

void main()
{
        double_complex z;
        int m,n;
        for (m = 0; m <= M; m++) 
        {
                // right-side streamline 
                for (n = -N; n <= N; n++) 
                {
                        z = g(m * dr + I * h(n * ds));
                        print_parts(z);
                }               
                // left-side streamline 
                for (n = -N; n <= N; n++) 
                {
                        z = g(-(m * dr) + I * h(n * ds));
                        print_parts(z);
                }               
        }
}

// To print  real(z)  and  imag(z)  to two files created with unique names ...

void print_parts(double_complex z)
{
        static FILE* restream;
        static FILE* imstream;
        static int isopen = 0;
        if (!isopen) 
        {
                char refile[L_tmpnam];
                char imfile[L_tmpnam];
                tmpnam(refile);
                tmpnam(imfile);
                restream = fopen(refile, "w");
                imstream = fopen(imfile, "w");
                printf("refile = %s\n", refile);
                printf("imfile = %s\n", imfile);
                isopen = 1;
        }
        fprintf(restream, " %17.10e", real(z));
        fprintf(imstream, " %17.10e", imag(z));
}