Subversion Repositories Nec2c

/*** Translated to the C language by N. Kyriazis  20 Aug 2003 ***

  Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,

  tape15,tape16,tape20,tape21)

  Numerical Electromagnetics Code (NEC2)  developed at Lawrence

  Livermore lab., Livermore, CA.  (contact G. Burke at 415-422-8414

  for problems with the NEC code. For problems with the vax implem-

  entation, contact J. Breakall at 415-422-8196 or E. Domning at 415

  422-5936)

  file created 4/11/80.

                                ***********Notice**********

 This computer code material was prepared as an account of work

 sponsored by the United States government.  Neither the United

 States nor the United States Department Of Energy, nor any of

 their employees, nor any of their contractors, subcontractors,

 or their employees, makes any warranty, express or implied, or

 assumes any legal liability or responsibility for the accuracy,

 completeness or usefulness of any information, apparatus, product

 or process disclosed, or represents that its use would not infringe

 privately-owned rights.

 ******************************************************************/

#include "nec2c.h"

/* common  /data/ */

extern data_t data;

/* common  /gnd/ */

extern gnd_t gnd;

/* common  /crnt/ */

extern crnt_t crnt;

/* common  /gwav/ */

extern gwav_t gwav;

/* common  /fpat/ */

extern fpat_t fpat;

/* pointers to input/output files */

extern FILE *input_fp, *output_fp, *plot_fp;

/* common  /save/ */

extern save_t save;

/* common  /plot/ */

extern plot_t plot;

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

/* ffld calculates the far zone radiated electric fields, */

/* the factor exp(j*k*r)/(r/lamda) not included */

void ffld (double thet, double phi, complex double *eth, complex double *eph)

{

        int k, i, ip, jump;

        double phx, phy, roz, rozs, thx, thy, thz, rox, roy;

        double tthet = 0., darg = 0., omega, el, sill, top, bot, a;

        double too, boo, b, c, d, rr, ri, arg, dr, rfl, rrz;

        complex double cix = CPLX_00, ciy = CPLX_00, ciz = CPLX_00;

        complex double exa, ccx = CPLX_00, ccy = CPLX_00, ccz = CPLX_00, cdp;

        complex double zrsin, rrv = CPLX_00, rrh = CPLX_00, rrv1 = CPLX_00;

        complex double rrh1 = CPLX_00, rrv2 = CPLX_00, rrh2 = CPLX_00;

        complex double tix, tiy, tiz, zscrn, ex = CPLX_00, ey = CPLX_00, ez = CPLX_00, gx, gy,

                                             gz;

        phx = -sinl (phi);

        phy = cosl (phi);

        roz = cosl (thet);

        rozs = roz;

        thx = roz * phy;

        thy = -roz * phx;

        thz = -sinl (thet);

        rox = -thz * phy;

        roy = thz * phx;

        jump = FALSE;

        if (data.n != 0)

        {

                /* loop for structure image if any */

                /* calculation of reflection coeffecients */

                for (k = 0; k < gnd.ksymp; k++)

                {

                        if (k != 0)

                        {

                                /* for perfect ground */

                                if (gnd.iperf == 1)

                                {

                                        rrv = -CPLX_10;

                                        rrh = -CPLX_10;

                                }

                                else

                                {

                                        /* for infinite planar ground */

                                        zrsin =

                                            csqrtl (1. - gnd.zrati * gnd.zrati * thz * thz);

                                        rrv = -(roz - gnd.zrati * zrsin) /

                                              (roz + gnd.zrati * zrsin);

                                        rrh = (gnd.zrati * roz - zrsin) /

                                              (gnd.zrati * roz + zrsin);

                                } /* if( gnd.iperf == 1) */

                                /* for the cliff problem, two reflction coefficients calculated

                                 */

                                if (gnd.ifar > 1)

                                {

                                        rrv1 = rrv;

                                        rrh1 = rrh;

                                        tthet = tanl (thet);

                                        if (gnd.ifar != 4)

                                        {

                                                zrsin = csqrtl (

                                                    1. - gnd.zrati2 * gnd.zrati2 * thz * thz);

                                                rrv2 = -(roz - gnd.zrati2 * zrsin) /

                                                       (roz + gnd.zrati2 * zrsin);

                                                rrh2 = (gnd.zrati2 * roz - zrsin) /

                                                       (gnd.zrati2 * roz + zrsin);

                                                darg = -TP * 2. * gnd.ch * roz;

                                        }

                                } /* if( gnd.ifar > 1) */

                                roz = -roz;

                                ccx = cix;

                                ccy = ciy;

                                ccz = ciz;

                        } /* if( k != 0 ) */

                        cix = CPLX_00;

                        ciy = CPLX_00;

                        ciz = CPLX_00;

                        /* loop over structure segments */

                        for (i = 0; i < data.n; i++)

                        {

                                omega =

                                    -(rox * data.cab[i] + roy * data.sab[i] +

                                      roz * data.salp[i]);

                                el = PI * data.si[i];

                                sill = omega * el;

                                top = el + sill;

                                bot = el - sill;

                                if (fabsl (omega) >= 1.0e-7)

                                        a = 2. * sinl (sill) / omega;

                                else

                                        a = (2. - omega * omega * el * el / 3.) * el;

                                if (fabsl (top) >= 1.0e-7)

                                        too = sinl (top) / top;

                                else

                                        too = 1. - top * top / 6.;

                                if (fabsl (bot) >= 1.0e-7)

                                        boo = sinl (bot) / bot;

                                else

                                        boo = 1. - bot * bot / 6.;

                                b = el * (boo - too);

                                c = el * (boo + too);

                                rr = a * crnt.air[i] + b * crnt.bii[i] + c * crnt.cir[i];

                                ri = a * crnt.aii[i] - b * crnt.bir[i] + c * crnt.cii[i];

                                arg =

                                    TP * (data.x[i] * rox + data.y[i] * roy + data.z[i] * roz);

                                if ((k != 1) || (gnd.ifar < 2))

                                {

                                        /* summation for far field integral */

                                        exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri);

                                        cix = cix + exa * data.cab[i];

                                        ciy = ciy + exa * data.sab[i];

                                        ciz = ciz + exa * data.salp[i];

                                        continue;

                                }

                                /* calculation of image contribution */

                                /* in cliff and ground screen problems */

                                /* specular point distance */

                                dr = data.z[i] * tthet;

                                d = dr * phy + data.x[i];

                                if (gnd.ifar == 2)

                                {

                                        if ((gnd.cl - d) > 0.)

                                        {

                                                rrv = rrv1;

                                                rrh = rrh1;

                                        }

                                        else

                                        {

                                                rrv = rrv2;

                                                rrh = rrh2;

                                                arg = arg + darg;

                                        }

                                } /* if( gnd.ifar == 2) */

                                else

                                {

                                        d = sqrtl (

                                            d * d +

                                            (data.y[i] - dr * phx) * (data.y[i] - dr * phx));

                                        if (gnd.ifar == 3)

                                        {

                                                if ((gnd.cl - d) > 0.)

                                                {

                                                        rrv = rrv1;

                                                        rrh = rrh1;

                                                }

                                                else

                                                {

                                                        rrv = rrv2;

                                                        rrh = rrh2;

                                                        arg = arg + darg;

                                                }

                                        } /* if( gnd.ifar == 3) */

                                        else

                                        {

                                                if ((gnd.scrwl - d) >= 0.)

                                                {

                                                        /* radial wire ground screen reflection

                                                         * coefficient */

                                                        d = d + gnd.t2;

                                                        zscrn = gnd.t1 * d * logl (d / gnd.t2);

                                                        zscrn = (zscrn * gnd.zrati) /

                                                                (ETA * gnd.zrati + zscrn);

                                                        zrsin = csqrtl (

                                                            1. - zscrn * zscrn * thz * thz);

                                                        rrv = (roz + zscrn * zrsin) /

                                                              (-roz + zscrn * zrsin);

                                                        rrh = (zscrn * roz + zrsin) /

                                                              (zscrn * roz - zrsin);

                                                } /* if(( gnd.scrwl- d) < 0.) */

                                                else

                                                {

                                                        if (gnd.ifar == 4)

                                                        {

                                                                rrv = rrv1;

                                                                rrh = rrh1;

                                                        } /* if( gnd.ifar == 4) */

                                                        else

                                                        {

                                                                if (gnd.ifar == 5)

                                                                        d = dr * phy +

                                                                            data.x[i];

                                                                if ((gnd.cl - d) > 0.)

                                                                {

                                                                        rrv = rrv1;

                                                                        rrh = rrh1;

                                                                }

                                                                else

                                                                {

                                                                        rrv = rrv2;

                                                                        rrh = rrh2;

                                                                        arg = arg + darg;

                                                                } /* if(( gnd.cl- d) > 0.) */

                                                        } /* if( gnd.ifar == 4) */

                                                } /* if(( gnd.scrwl- d) < 0.) */

                                        } /* if( gnd.ifar == 3) */

                                } /* if( gnd.ifar == 2) */

                                /* contribution of each image segment modified by */

                                /* reflection coef, for cliff and ground screen problems */

                                exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri);

                                tix = exa * data.cab[i];

                                tiy = exa * data.sab[i];

                                tiz = exa * data.salp[i];

                                cdp = (tix * phx + tiy * phy) * (rrh - rrv);

                                cix = cix + tix * rrv + cdp * phx;

                                ciy = ciy + tiy * rrv + cdp * phy;

                                ciz = ciz - tiz * rrv;

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

                        if (k == 0)

                                continue;

                        /* calculation of contribution of structure image for infinite ground

                         */

                        if (gnd.ifar < 2)

                        {

                                cdp = (cix * phx + ciy * phy) * (rrh - rrv);

                                cix = ccx + cix * rrv + cdp * phx;

                                ciy = ccy + ciy * rrv + cdp * phy;

                                ciz = ccz - ciz * rrv;

                        }

                        else

                        {

                                cix = cix + ccx;

                                ciy = ciy + ccy;

                                ciz = ciz + ccz;

                        }

                } /* for( k=0; k < gnd.ksymp; k++ ) */

                if (data.m > 0)

                        jump = TRUE;

                else

                {

                        *eth = (cix * thx + ciy * thy + ciz * thz) * CONST3;

                        *eph = (cix * phx + ciy * phy) * CONST3;

                        return;

                }

        } /* if( n != 0) */

        if (!jump)

        {

                cix = CPLX_00;

                ciy = CPLX_00;

                ciz = CPLX_00;

        }

        /* electric field components */

        roz = rozs;

        rfl = -1.;

        for (ip = 0; ip < gnd.ksymp; ip++)

        {

                rfl = -rfl;

                rrz = roz * rfl;

                fflds (rox, roy, rrz, &crnt.cur[data.n], &gx, &gy, &gz);

                if (ip != 1)

                {

                        ex = gx;

                        ey = gy;

                        ez = gz;

                        continue;

                }

                if (gnd.iperf == 1)

                {

                        gx = -gx;

                        gy = -gy;

                        gz = -gz;

                }

                else

                {

                        rrv = csqrtl (1. - gnd.zrati * gnd.zrati * thz * thz);

                        rrh = gnd.zrati * roz;

                        rrh = (rrh - rrv) / (rrh + rrv);

                        rrv = gnd.zrati * rrv;

                        rrv = -(roz - rrv) / (roz + rrv);

                        *eth = (gx * phx + gy * phy) * (rrh - rrv);

                        gx = gx * rrv + *eth * phx;

                        gy = gy * rrv + *eth * phy;

                        gz = gz * rrv;

                } /* if( gnd.iperf == 1) */

                ex = ex + gx;

                ey = ey + gy;

                ez = ez - gz;

        } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */

        ex = ex + cix * CONST3;

        ey = ey + ciy * CONST3;

        ez = ez + ciz * CONST3;

        *eth = ex * thx + ey * thy + ez * thz;

        *eph = ex * phx + ey * phy;

        return;

}

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

/* calculates the xyz components of the electric */

/* field due to surface currents */

void fflds (

    double rox,

    double roy,

    double roz,

    complex double *scur,

    complex double *ex,

    complex double *ey,

    complex double *ez)

{

        double *xs, *ys, *zs, *s;

        int j, i, k;

        double arg;

        complex double ct;

        xs = data.px;

        ys = data.py;

        zs = data.pz;

        s = data.pbi;

        *ex = CPLX_00;

        *ey = CPLX_00;

        *ez = CPLX_00;

        i = -1;

        for (j = 0; j < data.m; j++)

        {

                i++;

                arg = TP * (rox * xs[i] + roy * ys[i] + roz * zs[i]);

                ct = cmplx (cosl (arg) * s[i], sinl (arg) * s[i]);

                k = 3 * j;

                *ex += scur[k] * ct;

                *ey += scur[k + 1] * ct;

                *ez += scur[k + 2] * ct;

        }

        ct = rox * *ex + roy * *ey + roz * *ez;

        *ex = CONST4 * (ct * rox - *ex);

        *ey = CONST4 * (ct * roy - *ey);

        *ez = CONST4 * (ct * roz - *ez);

        return;

}

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

/* gfld computes the radiated field including ground wave. */

void gfld (

    double rho,

    double phi,

    double rz,

    complex double *eth,

    complex double *epi,

    complex double *erd,

    complex double ux,

    int ksymp)

{

        int i, k;

        double b, r, thet, arg, phx, phy, rx, ry, dx, dy, dz, rix, riy, rhs, rhp;

        double rhx, rhy, calp, cbet, sbet, cph, sph, el, rfl, riz, thx, thy, thz;

        double rxyz, rnx, rny, rnz, omega, sill, top, bot, a, too, boo, c, rr, ri;

        complex double cix, ciy, ciz, exa, erv;

        complex double ezv, erh, eph, ezh, ex, ey;

        r = sqrtl (rho * rho + rz * rz);

        if ((ksymp == 1) || (cabs (ux) > .5) || (r > 1.e5))

        {

                /* computation of space wave only */

                if (rz >= 1.0e-20)

                        thet = atanl (rho / rz);

                else

                        thet = PI * .5;

                ffld (thet, phi, eth, epi);

                arg = -TP * r;

                exa = cmplx (cosl (arg), sinl (arg)) / r;

                *eth = *eth * exa;

                *epi = *epi * exa;

                *erd = CPLX_00;

                return;

        } /* if( (ksymp == 1) && (cabs(ux) > .5) && (r > 1.e5) ) */

        /* computation of space and ground waves. */

        gwav.u = ux;

        gwav.u2 = gwav.u * gwav.u;

        phx = -sinl (phi);

        phy = cosl (phi);

        rx = rho * phy;

        ry = -rho * phx;

        cix = CPLX_00;

        ciy = CPLX_00;

        ciz = CPLX_00;

        /* summation of field from individual segments */

        for (i = 0; i < data.n; i++)

        {

                dx = data.cab[i];

                dy = data.sab[i];

                dz = data.salp[i];

                rix = rx - data.x[i];

                riy = ry - data.y[i];

                rhs = rix * rix + riy * riy;

                rhp = sqrtl (rhs);

                if (rhp >= 1.0e-6)

                {

                        rhx = rix / rhp;

                        rhy = riy / rhp;

                }

                else

                {

                        rhx = 1.;

                        rhy = 0.;

                }

                calp = 1. - dz * dz;

                if (calp >= 1.0e-6)

                {

                        calp = sqrtl (calp);

                        cbet = dx / calp;

                        sbet = dy / calp;

                        cph = rhx * cbet + rhy * sbet;

                        sph = rhy * cbet - rhx * sbet;

                }

                else

                {

                        cph = rhx;

                        sph = rhy;

                }

                el = PI * data.si[i];

                rfl = -1.;

                /* integration of (current)*(phase factor) over segment and image for */

                /* constant, sine, and cosine current distributions */

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

                {

                        rfl = -rfl;

                        riz = rz - data.z[i] * rfl;

                        rxyz = sqrtl (rix * rix + riy * riy + riz * riz);

                        rnx = rix / rxyz;

                        rny = riy / rxyz;

                        rnz = riz / rxyz;

                        omega = -(rnx * dx + rny * dy + rnz * dz * rfl);

                        sill = omega * el;

                        top = el + sill;

                        bot = el - sill;

                        if (fabsl (omega) >= 1.0e-7)

                                a = 2. * sinl (sill) / omega;

                        else

                                a = (2. - omega * omega * el * el / 3.) * el;

                        if (fabsl (top) >= 1.0e-7)

                                too = sinl (top) / top;

                        else

                                too = 1. - top * top / 6.;

                        if (fabsl (bot) >= 1.0e-7)

                                boo = sinl (bot) / bot;

                        else

                                boo = 1. - bot * bot / 6.;

                        b = el * (boo - too);

                        c = el * (boo + too);

                        rr = a * crnt.air[i] + b * crnt.bii[i] + c * crnt.cir[i];

                        ri = a * crnt.aii[i] - b * crnt.bir[i] + c * crnt.cii[i];

                        arg = TP * (data.x[i] * rnx + data.y[i] * rny + data.z[i] * rnz * rfl);

                        exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri) / TP;

                        if (k != 1)

                        {

                                gwav.xx1 = exa;

                                gwav.r1 = rxyz;

                                gwav.zmh = riz;

                                continue;

                        }

                        gwav.xx2 = exa;

                        gwav.r2 = rxyz;

                        gwav.zph = riz;

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

                /* call subroutine to compute the field */

                /* of segment including ground wave. */

                gwave (&erv, &ezv, &erh, &ezh, &eph);

                erh = erh * cph * calp + erv * dz;

                eph = eph * sph * calp;

                ezh = ezh * cph * calp + ezv * dz;

                ex = erh * rhx - eph * rhy;

                ey = erh * rhy + eph * rhx;

                cix = cix + ex;

                ciy = ciy + ey;

                ciz = ciz + ezh;

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

        arg = -TP * r;

        exa = cmplx (cosl (arg), sinl (arg));

        cix = cix * exa;

        ciy = ciy * exa;

        ciz = ciz * exa;

        rnx = rx / r;

        rny = ry / r;

        rnz = rz / r;

        thx = rnz * phy;

        thy = -rnz * phx;

        thz = -rho / r;

        *eth = cix * thx + ciy * thy + ciz * thz;

        *epi = cix * phx + ciy * phy;

        *erd = cix * rnx + ciy * rny + ciz * rnz;

        return;

}

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

/* compute radiation pattern, gain, normalized gain */

void rdpat (void)

{

        char *hpol[3] = {"LINEAR", "RIGHT ", "LEFT  "};

        char *igtp[2] = {"- POWER GAINS -", "- DIRECTIVE GAINS -"};

        char *igax[4] = {" MAJOR", " MINOR", "VERT. ", "HOR.  "};

        char *igntp[5] = {

            " MAJOR AXIS", "  MINOR AXIS", "    VERTICAL", "  HORIZONTAL", "       TOTAL "};

        char *hclif = NULL, *isens;

        int i, j, jump, itmp1, itmp2, kth, kph, itmp3, itmp4;

        double exrm = 0., exra = 0., prad, gcon, gcop, gmax, pint, tmp1, tmp2;

        double phi, pha, thet, tha, erdm = 0., erda = 0., ethm2, ethm, *gain = NULL;

        double etha, ephm2, ephm, epha, tilta, emajr2, eminr2, axrat;

        double dfaz, dfaz2, cdfaz, tstor1 = 0., tstor2, stilta, gnmj;

        double gnmn, gnv, gnh, gtot, tmp3, tmp4, da, tmp5, tmp6;

        complex double eth, eph, erd;

        /* Allocate memory to gain buffer */

        if (fpat.inor > 0)

                mem_alloc ((void *) &gain, fpat.nth * fpat.nph * sizeof (double));

        if (gnd.ifar > 1)

        {

                fprintf (

                    output_fp,

                    "\n\n\n\n"

                    "                               "

                    "- - - FAR FIELD GROUND PARAMETERS - - -\n");

                jump = FALSE;

                if (gnd.ifar > 3)

                {

                        fprintf (

                            output_fp,

                            "\n"

                            "                               "

                            "--- RADIAL WIRE GROUND SCREEN ---\n"

                            "                               "

                            "NUM OF WIRES= %d\n"

                            "                               "

                            "WIRE LENGTH= %8.2f METERS\n"

                            "                               "

                            "WIRE RADIUS= %10.3E METERS",

                            gnd.nradl,

                            save.scrwlt,

                            save.scrwrt);

                        if (gnd.ifar == 4)

                                jump = TRUE;

                } /* if( gnd.ifar > 3) */

                if (!jump)

                {

                        if ((gnd.ifar == 2) || (gnd.ifar == 5))

                                hclif = "LINEAR";

                        if ((gnd.ifar == 3) || (gnd.ifar == 6))

                                hclif = "CIRCULAR";

                        gnd.cl = fpat.clt / data.wlam;

                        gnd.ch = fpat.cht / data.wlam;

                        gnd.zrati2 =

                            csqrtl (1. / cmplx (fpat.epsr2, -fpat.sig2 * data.wlam * 59.96));

                        fprintf (

                            output_fp,

                            "\n"

                            "                               "

                            "- - %s CLIFF - -\n"

                            "                               "

                            "EDGE DISTANCE= %9.2f METERS\n"

                            "                               "

                            "       HEIGHT= %9.2f METERS\n"

                            "                               "

                            "--- SECOND MEDIUM ---\n"

                            "                               "

                            "RELATIVE DIELECTRIC CONST= %10.3f\n"

                            "                               "

                            "      GROUND CONDUCTIVITY= %10.3f MHOS",

                            hclif,

                            fpat.clt,

                            fpat.cht,

                            fpat.epsr2,

                            fpat.sig2);

                } /* if( ! jump ) */

        } /* if( gnd.ifar > 1) */

        if (gnd.ifar == 1)

        {

                fprintf (

                    output_fp,

                    "\n\n\n"

                    "                             "

                    "------- RADIATED FIELDS NEAR GROUND --------\n\n"

                    "    ------- LOCATION -------     --- E(THETA) ---    "

                    " ---- E(PHI) ----    --- E(RADIAL) ---\n"

                    "      RHO    PHI        Z           MAG    PHASE     "

                    "    MAG    PHASE        MAG     PHASE \n"

                    "    METERS   DEG.     METERS      VOLTS/M DEGREES   "

                    "   VOLTS/M   DEG.     VOLTS/M    DEG. ");

        }

        else

        {

                itmp1 = 2 * fpat.iax;

                itmp2 = itmp1 + 1;

                fprintf (

                    output_fp,

                    "\n\n\n\n"

                    "                                                "

                    "- - - RADIATION PATTERNS - - -\n");

                if (fpat.rfld >= 1.0e-20)

                {

                        exrm = 1. / fpat.rfld;

                        exra = fpat.rfld / data.wlam;

                        exra = -360. * (exra - floorl (exra));

                        fprintf (

                            output_fp,

                            "\n"

                            "                             "

                            "RANGE: %13.6E METERS\n"

                            "                             "

                            "EXP(-JKR)/R: %12.5E AT PHASE: %7.2f DEGREES\n",

                            fpat.rfld,

                            exrm,

                            exra);

                }

                fprintf (

                    output_fp,

                    "\n"

                    "  - - ANGLES - -   %23s       - - - POLARIZATION - - -  "

                    "  - - - E(THETA) - - -    - - - E(PHI) - - -\n"

                    "  THETA     PHI        %6s  %6s  TOTAL      AXIAL     TILT  "

                    " SENSE     MAGNITUDE    PHASE      MAGNITUDE    PHASE \n"

                    " DEGREES  DEGREES       DB      DB      DB        RATIO     "

                    "DEG.              VOLTS/M    DEGREES      VOLTS/M    DEGREES",

                    igtp[fpat.ipd],

                    igax[itmp1],

                    igax[itmp2]);

        } /* if( gnd.ifar == 1) */

        if ((fpat.ixtyp == 0) || (fpat.ixtyp == 5))

        {

                gcop = data.wlam * data.wlam * 2. * PI / (376.73 * fpat.pinr);

                prad = fpat.pinr - fpat.ploss - fpat.pnlr;

                gcon = gcop;

                if (fpat.ipd != 0)

                        gcon = gcon * fpat.pinr / prad;

        }

        else if (fpat.ixtyp == 4)

        {

                fpat.pinr = 394.51 * fpat.xpr6 * fpat.xpr6 * data.wlam * data.wlam;

                gcop = data.wlam * data.wlam * 2. * PI / (376.73 * fpat.pinr);

                prad = fpat.pinr - fpat.ploss - fpat.pnlr;

                gcon = gcop;