Subversion Repositories Nec2c

/* last change:  pgm   8 nov 2000    1:04 pm

 program somnec(input,output,tape21)

 program to generate nec interpolation grids for fields due to

 ground.  field components are computed by numerical evaluation

 of modified sommerfeld integrals.

 somnec2d is a double precision version of somnec for use with

 nec2d.  an alternate version (somnec2sd) is also provided in which

 computation is in single precision but the output file is written

 in double precision for use with nec2d.  somnec2sd runs about twic

 as fast as the full double precision somnec2d.  the difference

 between nec2d results using a for021 file from this code rather

 than from somnec2sd was insignficant in the cases tested.

 changes made by j bergervoet, 31-5-95:

 parameter 0. --> 0.d0 in calling of routine test

 status of output files set to 'unknown' */

#include "nec2c.h"

/* common /evlcom/ */

static int jh;

static double ck2, ck2sq, tkmag, tsmag, ck1r, zph, rho;

static complex double ct1, ct2, ct3, ck1, ck1sq, cksm;

/* common /cntour/ */

static complex double a, b;

/*common  /ggrid/ */

ggrid_t ggrid;

/*-----------------------------------------------------------------------*/

/* This is the "main" of somnec */

void somnec (double epr, double sig, double fmhz)

{

        int k, nth, ith, irs, ir, nr;

        double tim, wlam, tst, dr, dth, r, rk, thet, tfac1, tfac2;

        complex double erv, ezv, erh, eph, cl1, cl2, con;

        if (sig >= 0.)

        {

                wlam = CVEL / fmhz;

                ggrid.epscf = cmplx (epr, -sig * wlam * 59.96);

        }

        else

                ggrid.epscf = cmplx (epr, sig);

        secnds (&tst);

        ck2 = TP;

        ck2sq = ck2 * ck2;

        /* sommerfeld integral evaluation uses exp(-jwt), nec uses exp(+jwt), */

        /* hence need conjg(ggrid.epscf).  conjugate of fields occurs in subroutine */

        /* evlua. */

        ck1sq = ck2sq * conj (ggrid.epscf);

        ck1 = csqrtl (ck1sq);

        ck1r = creal (ck1);

        tkmag = 100. * cabs (ck1);

        tsmag = 100. * ck1 * conj (ck1);

        cksm = ck2sq / (ck1sq + ck2sq);

        ct1 = .5 * (ck1sq - ck2sq);

        erv = ck1sq * ck1sq;

        ezv = ck2sq * ck2sq;

        ct2 = .125 * (erv - ezv);

        erv *= ck1sq;

        ezv *= ck2sq;

        ct3 = .0625 * (erv - ezv);

        /* loop over 3 grid regions */

        for (k = 0; k < 3; k++)

        {

                nr = ggrid.nxa[k];

                nth = ggrid.nya[k];

                dr = ggrid.dxa[k];

                dth = ggrid.dya[k];

                r = ggrid.xsa[k] - dr;

                irs = 1;

                if (k == 0)

                {

                        r = ggrid.xsa[k];

                        irs = 2;

                }

                /*  loop over r.  (r=sqrtl(rho**2 + (z+h)**2)) */

                for (ir = irs - 1; ir < nr; ir++)

                {

                        r += dr;

                        thet = ggrid.ysa[k] - dth;

                        /* loop over theta.  (theta=atan((z+h)/rho)) */

                        for (ith = 0; ith < nth; ith++)

                        {

                                thet += dth;

                                rho = r * cosl (thet);

                                zph = r * sinl (thet);

                                if (rho < 1.e-7)

                                        rho = 1.e-8;

                                if (zph < 1.e-7)

                                        zph = 0.;

                                evlua (&erv, &ezv, &erh, &eph);

                                rk = ck2 * r;

                                con = -CONST1 * r / cmplx (cosl (rk), -sinl (rk));

                                switch (k)

                                {

                                case 0:

                                        ggrid.ar1[ir + ith * 11 + 0] = erv * con;

                                        ggrid.ar1[ir + ith * 11 + 110] = ezv * con;

                                        ggrid.ar1[ir + ith * 11 + 220] = erh * con;

                                        ggrid.ar1[ir + ith * 11 + 330] = eph * con;

                                        break;

                                case 1:

                                        ggrid.ar2[ir + ith * 17 + 0] = erv * con;

                                        ggrid.ar2[ir + ith * 17 + 85] = ezv * con;

                                        ggrid.ar2[ir + ith * 17 + 170] = erh * con;

                                        ggrid.ar2[ir + ith * 17 + 255] = eph * con;

                                        break;

                                case 2:

                                        ggrid.ar3[ir + ith * 9 + 0] = erv * con;

                                        ggrid.ar3[ir + ith * 9 + 72] = ezv * con;

                                        ggrid.ar3[ir + ith * 9 + 144] = erh * con;

                                        ggrid.ar3[ir + ith * 9 + 216] = eph * con;

                                } /* switch( k ) */

                        } /* for( ith = 0; ith < nth; ith++ ) */

                } /* for( ir = irs-1; ir < nr; ir++; ) */

        } /* for( k = 0; k < 3; k++; ) */

        /* fill grid 1 for r equal to zero. */

        cl2 = -CONST4 * (ggrid.epscf - 1.) / (ggrid.epscf + 1.);

        cl1 = cl2 / (ggrid.epscf + 1.);

        ezv = ggrid.epscf * cl1;

        thet = -dth;

        nth = ggrid.nya[0];

        for (ith = 0; ith < nth; ith++)

        {

                thet += dth;

                if ((ith + 1) != nth)

                {

                        tfac2 = cosl (thet);

                        tfac1 = (1. - sinl (thet)) / tfac2;

                        tfac2 = tfac1 / tfac2;

                        erv = ggrid.epscf * cl1 * tfac1;

                        erh = cl1 * (tfac2 - 1.) + cl2;

                        eph = cl1 * tfac2 - cl2;

                }

                else

                {

                        erv = 0.;

                        erh = cl2 - .5 * cl1;

                        eph = -erh;

                }

                ggrid.ar1[0 + ith * 11 + 0] = erv;

                ggrid.ar1[0 + ith * 11 + 110] = ezv;

                ggrid.ar1[0 + ith * 11 + 220] = erh;

                ggrid.ar1[0 + ith * 11 + 330] = eph;

        }

        secnds (&tim);

        tim -= tst;

        return;

}

/*-----------------------------------------------------------------------*/

/* bessel evaluates the zero-order bessel function */

/* and its derivative for complex argument z. */

void bessel (complex double z, complex double *j0, complex double *j0p)

{

        int k, i, ib, iz, miz;

        static int m[101], init = FALSE;

        static double a1[25], a2[25];

        double tst, zms;

        complex double p0z, p1z, q0z, q1z, zi, zi2, zk, cz, sz, j0x = CPLX_00, j0px = CPLX_00;

        /* initialization of constants */

        if (!init)

        {

                for (k = 1; k <= 25; k++)

                {

                        i = k - 1;

                        a1[i] = -.25 / (k * k);

                        a2[i] = 1.0 / (k + 1.0);

                }

                for (i = 1; i <= 101; i++)

                {

                        tst = 1.0;

                        for (k = 0; k < 24; k++)

                        {

                                init = k;

                                tst *= -i * a1[k];

                                if (tst < 1.0e-6)

                                        break;

                        }

                        m[i - 1] = init + 1;

                } /* for( i = 1; i<= 101; i++ ) */

                init = TRUE;

        } /* if(init == 0) */

        zms = z * conj (z);

        if (zms <= 1.e-12)

        {

                *j0 = CPLX_10;

                *j0p = -.5 * z;

                return;

        }

        ib = 0;

        if (zms <= 37.21)

        {

                if (zms > 36.)

                        ib = 1;

                /* series expansion */

                iz = zms;

                miz = m[iz];

                *j0 = CPLX_10;

                *j0p = *j0;

                zk = *j0;

                zi = z * z;

                for (k = 0; k < miz; k++)

                {

                        zk *= a1[k] * zi;

                        *j0 += zk;

                        *j0p += a2[k] * zk;

                }

                *j0p *= -.5 * z;

                if (ib == 0)

                        return;

                j0x = *j0;

                j0px = *j0p;

        }

        /* asymptotic expansion */

        zi = 1. / z;

        zi2 = zi * zi;

        p0z = 1. + (P20 * zi2 - P10) * zi2;

        p1z = 1. + (P11 - P21 * zi2) * zi2;

        q0z = (Q20 * zi2 - Q10) * zi;

        q1z = (Q11 - Q21 * zi2) * zi;

        zk = cexp (CPLX_01 * (z - POF));

        zi2 = 1. / zk;

        cz = .5 * (zk + zi2);

        sz = CPLX_01 * .5 * (zi2 - zk);

        zk = C3 * csqrtl (zi);

        *j0 = zk * (p0z * cz - q0z * sz);

        *j0p = -zk * (p1z * sz + q1z * cz);

        if (ib == 0)

                return;

        zms = cosl ((sqrtl (zms) - 6.) * PI10);

        *j0 = .5 * (j0x * (1. + zms) + *j0 * (1. - zms));

        *j0p = .5 * (j0px * (1. + zms) + *j0p * (1. - zms));

        return;

}

/*-----------------------------------------------------------------------*/

/* evlua controls the integration contour in the complex */

/* lambda plane for evaluation of the sommerfeld integrals */

void evlua (complex double *erv, complex double *ezv, complex double *erh, complex double *eph)

{

        int i, jump;

        static double del, slope, rmis;

        static complex double cp1, cp2, cp3, bk, delta, delta2, sum[6], ans[6];

        del = zph;

        if (rho > del)

                del = rho;

        if (zph >= 2. * rho)

        {

                /* bessel function form of sommerfeld integrals */

                jh = 0;

                a = CPLX_00;

                del = 1. / del;

                if (del > tkmag)

                {

                        b = cmplx (.1 * tkmag, -.1 * tkmag);

                        rom1 (6, sum, 2);

                        a = b;

                        b = cmplx (del, -del);

                        rom1 (6, ans, 2);

                        for (i = 0; i < 6; i++)

                                sum[i] += ans[i];

                }

                else

                {

                        b = cmplx (del, -del);

                        rom1 (6, sum, 2);

                }

                delta = PTP * del;

                gshank (b, delta, ans, 6, sum, 0, b, b);

                ans[5] *= ck1;

                /* conjugate since nec uses exp(+jwt) */

                *erv = conj (ck1sq * ans[2]);

                *ezv = conj (ck1sq * (ans[1] + ck2sq * ans[4]));

                *erh = conj (ck2sq * (ans[0] + ans[5]));

                *eph = -conj (ck2sq * (ans[3] + ans[5]));

                return;

        } /* if(zph >= 2.*rho) */

        /* hankel function form of sommerfeld integrals */

        jh = 1;

        cp1 = cmplx (0.0, .4 * ck2);

        cp2 = cmplx (.6 * ck2, -.2 * ck2);

        cp3 = cmplx (1.02 * ck2, -.2 * ck2);

        a = cp1;

        b = cp2;

        rom1 (6, sum, 2);

        a = cp2;

        b = cp3;

        rom1 (6, ans, 2);

        for (i = 0; i < 6; i++)

                sum[i] = -(sum[i] + ans[i]);

        /* path from imaginary axis to -infinity */

        if (zph > .001 * rho)

                slope = rho / zph;

        else

                slope = 1000.;

        del = PTP / del;

        delta = cmplx (-1.0, slope) * del / sqrtl (1. + slope * slope);

        delta2 = -conj (delta);

        gshank (cp1, delta, ans, 6, sum, 0, bk, bk);

        rmis = rho * (creal (ck1) - ck2);

        jump = FALSE;

        if ((rmis >= 2. * ck2) && (rho >= 1.e-10))

        {

                if (zph >= 1.e-10)

                {

                        bk = cmplx (-zph, rho) * (ck1 - cp3);

                        rmis = -creal (bk) / fabsl (cimag (bk));

                        if (rmis > 4. * rho / zph)

                                jump = TRUE;

                }

                if (!jump)

                {

                        /* integrate up between branch cuts, then to + infinity */

                        cp1 = ck1 - (.1 + .2fj);

                        cp2 = cp1 + .2;

                        bk = cmplx (0., del);

                        gshank (cp1, bk, sum, 6, ans, 0, bk, bk);

                        a = cp1;

                        b = cp2;

                        rom1 (6, ans, 1);

                        for (i = 0; i < 6; i++)

                                ans[i] -= sum[i];

                        gshank (cp3, bk, sum, 6, ans, 0, bk, bk);

                        gshank (cp2, delta2, ans, 6, sum, 0, bk, bk);

                }

                jump = TRUE;

        } /* if( (rmis >= 2.*ck2) || (rho >= 1.e-10) ) */

        else

                jump = FALSE;

        if (!jump)

        {

                /* integrate below branch points, then to + infinity */

                for (i = 0; i < 6; i++)

                        sum[i] = -ans[i];

                rmis = creal (ck1) * 1.01;

                if ((ck2 + 1.) > rmis)

                        rmis = ck2 + 1.;

                bk = cmplx (rmis, .99 * cimag (ck1));

                delta = bk - cp3;

                delta *= del / cabs (delta);

                gshank (cp3, delta, ans, 6, sum, 1, bk, delta2);

        } /* if( ! jump ) */

        ans[5] *= ck1;

        /* conjugate since nec uses exp(+jwt) */

        *erv = conj (ck1sq * ans[2]);

        *ezv = conj (ck1sq * (ans[1] + ck2sq * ans[4]));

        *erh = conj (ck2sq * (ans[0] + ans[5]));

        *eph = -conj (ck2sq * (ans[3] + ans[5]));

        return;

}

/*-----------------------------------------------------------------------*/

/* fbar is sommerfeld attenuation function for numerical distance p */

void fbar (complex double p, complex double *fbar)

{

        int i, minus;

        double tms, sms;

        complex double z, zs, sum, pow, term;

        z = CPLX_01 * csqrtl (p);

        if (cabs (z) <= 3.)

        {

                /* series expansion */

                zs = z * z;

                sum = z;

                pow = z;

                for (i = 1; i <= 100; i++)

                {

                        pow = -pow * zs / (double) i;

                        term = pow / (2. * i + 1.);

                        sum = sum + term;

                        tms = creal (term * conj (term));

                        sms = creal (sum * conj (sum));

                        if (tms / sms < ACCS)

                                break;

                }

                *fbar = 1. - (1. - sum * TOSP) * z * cexp (zs) * SP;

        } /* if( cabs( z) <= 3.) */

        /* asymptotic expansion */

        if (creal (z) < 0.)

        {

                minus = 1;

                z = -z;

        }

        else

                minus = 0;

        zs = .5 / (z * z);

        sum = CPLX_00;

        term = CPLX_10;

        for (i = 1; i <= 6; i++)

        {

                term = -term * (2. * i - 1.) * zs;

                sum += term;

        }

        if (minus == 1)

                sum -= 2. * SP * z * cexp (z * z);

        *fbar = -sum;

}

/*-----------------------------------------------------------------------*/

/* gshank integrates the 6 sommerfeld integrals from start to */

/* infinity (until convergence) in lambda.  at the break point, bk, */

/* the step increment may be changed from dela to delb.  shank's */

/* algorithm to accelerate convergence of a slowly converging series */

/* is used */

void gshank (

    complex double start,

    complex double dela,

    complex double *sum,

    int nans,

    complex double *seed,

    int ibk,

    complex double bk,

    complex double delb)

{

        int ibx, j, i, jm, intx, inx, brk = 0;

        static double rbk, amg, den, denm;

        complex double a1, a2, as1, as2, del, aa;

        complex double q1[6][20], q2[6][20], ans1[6], ans2[6];

        rbk = creal (bk);

        del = dela;

        if (ibk == 0)

                ibx = 1;

        else

                ibx = 0;

        for (i = 0; i < nans; i++)

                ans2[i] = seed[i];

        b = start;

        for (intx = 1; intx <= MAXH; intx++)

        {

                inx = intx - 1;

                a = b;

                b += del;

                if ((ibx == 0) && (creal (b) >= rbk))

                {

                        /* hit break point.  reset seed and start over. */

                        ibx = 1;

                        b = bk;

                        del = delb;

                        rom1 (nans, sum, 2);

                        if (ibx != 2)

                        {

                                for (i = 0; i < nans; i++)

                                        ans2[i] += sum[i];

                                intx = 0;

                                continue;

                        }

                        for (i = 0; i < nans; i++)

                                ans2[i] = ans1[i] + sum[i];

                        intx = 0;

                        continue;

                } /* if( (ibx == 0) && (creal(b) >= rbk) ) */

                rom1 (nans, sum, 2);

                for (i = 0; i < nans; i++)

                        ans1[i] = ans2[i] + sum[i];

                a = b;

                b += del;

                if ((ibx == 0) && (creal (b) >= rbk))

                {

                        /* hit break point.  reset seed and start over. */

                        ibx = 2;

                        b = bk;

                        del = delb;

                        rom1 (nans, sum, 2);

                        if (ibx != 2)

                        {

                                for (i = 0; i < nans; i++)

                                        ans2[i] += sum[i];

                                intx = 0;

                                continue;

                        }

                        for (i = 0; i < nans; i++)

                                ans2[i] = ans1[i] + sum[i];

                        intx = 0;

                        continue;

                } /* if( (ibx == 0) && (creal(b) >= rbk) ) */

                rom1 (nans, sum, 2);

                for (i = 0; i < nans; i++)

                        ans2[i] = ans1[i] + sum[i];

                den = 0.;

                for (i = 0; i < nans; i++)

                {

                        as1 = ans1[i];

                        as2 = ans2[i];

                        if (intx >= 2)

                        {

                                for (j = 1; j < intx; j++)

                                {

                                        jm = j - 1;

                                        aa = q2[i][jm];

                                        a1 = q1[i][jm] + as1 - 2. * aa;

                                        if ((creal (a1) != 0.) || (cimag (a1) != 0.))

                                        {

                                                a2 = aa - q1[i][jm];

                                                a1 = q1[i][jm] - a2 * a2 / a1;

                                        }

                                        else

                                                a1 = q1[i][jm];

                                        a2 = aa + as2 - 2. * as1;

                                        if ((creal (a2) != 0.) || (cimag (a2) != 0.))

                                                a2 = aa - (as1 - aa) * (as1 - aa) / a2;

                                        else

                                                a2 = aa;

                                        q1[i][jm] = as1;

                                        q2[i][jm] = as2;

                                        as1 = a1;

                                        as2 = a2;

                                } /* for( j = 1; i < intx; i++ ) */

                        } /* if(intx >= 2) */

                        q1[i][intx - 1] = as1;

                        q2[i][intx - 1] = as2;

                        amg = fabsl (creal (as2)) + fabsl (cimag (as2));

                        if (amg > den)

                                den = amg;

                } /* for( i = 0; i < nans; i++ ) */

                denm = 1.e-3 * den * CRIT;

                jm = intx - 3;

                if (jm < 1)

                        jm = 1;

                for (j = jm - 1; j < intx; j++)

                {

                        brk = FALSE;

                        for (i = 0; i < nans; i++)

                        {

                                a1 = q2[i][j];

                                den = (fabsl (creal (a1)) + fabsl (cimag (a1))) * CRIT;

                                if (den < denm)

                                        den = denm;

                                a1 = q1[i][j] - a1;

                                amg = fabsl (creal (a1) + fabsl (cimag (a1)));

                                if (amg > den)

                                {

                                        brk = TRUE;

                                        break;

                                }

                        } /* for( i = 0; i < nans; i++ ) */

                        if (brk)

                                break;

                } /* for( j = jm-1; j < intx; j++ ) */

                if (!brk)

                {

                        for (i = 0; i < nans; i++)

                                sum[i] = .5 * (q1[i][inx] + q2[i][inx]);

                        return;

                }

        } /* for( intx = 1; intx <= maxh; intx++ ) */

        /* No convergence */

        abort_on_error (-6);

}

/*-----------------------------------------------------------------------*/

/* hankel evaluates hankel function of the first kind,   */

/* order zero, and its derivative for complex argument z */

void hankel (complex double z, complex double *h0, complex double *h0p)

{

        int i, k, ib, iz, miz;

        static int m[101], init = FALSE;

        static double a1[25], a2[25], a3[25], a4[25], psi, tst, zms;

        complex double clogz, j0, j0p, p0z, p1z, q0z, q1z, y0 = CPLX_00, y0p = CPLX_00, zi,

                                                           zi2, zk;

        /* initialization of constants */

        if (!init)

        {

                psi = -GAMMA;

                for (k = 1; k <= 25; k++)

                {

                        i = k - 1;

                        a1[i] = -.25 / (k * k);

                        a2[i] = 1.0 / (k + 1.0);

                        psi += 1.0 / k;

                        a3[i] = psi + psi;

                        a4[i] = (psi + psi + 1.0 / (k + 1.0)) / (k + 1.0);

                }

                for (i = 1; i <= 101; i++)

                {

                        tst = 1.0;

                        for (k = 0; k < 24; k++)

                        {

                                init = k;

                                tst *= -i * a1[k];

                                if (tst * a3[k] < 1.e-6)

                                        break;

                        }

                        m[i - 1] = init + 1;

                }

                init = TRUE;

        } /* if( ! init ) */

        zms = z * conj (z);

        if (zms == 0.)

                abort_on_error (-7);

        ib = 0;

        if (zms <= 16.81)

        {

                if (zms > 16.)

                        ib = 1;

                /* series expansion */

                iz = zms;

                miz = m[iz];

                j0 = CPLX_10;

                j0p = j0;

                y0 = CPLX_00;

                y0p = y0;

                zk = j0;