#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define RAD     (180./3.14159)

static int kmdl=-1;
static double amr, cx, cy, cz, cxy, cyz, czx;
static double xr[2], yr[2], zr[2], wr, hsv, dcx, dcy, dcz;

void magafd(double the, double phi, double dip, double dec)
{   double cthe, cdip, px, py, pz, ex, ey, ez;
    cthe = cos(the/RAD); cdip = cos(dip/RAD);
    px = cthe*cos(phi/RAD); py = cthe*sin(phi/RAD); pz = sin(the/RAD);
    ex = cdip*cos(dec/RAD); ey = cdip*sin(dec/RAD); ez = sin(dip/RAD);
    cx = ex*px; cy = ey*py; cz = ez*pz;
    cxy = (ex*py+ey*px)/2.; cyz = (ey*pz+ez*py)/2.; czx = (ez*px+ex*pz)/2.;
}

void mpoint(double am, double xc, double yc, double zc, double vol)
{   amr = am * 100.; hsv = vol; kmdl = 0;
    xr[0] = xc; yr[0] = yc; zr[0] = zc;
}

void mvline(double am, double xc, double yc, double za, double w, double area)
{   amr = am * 100.; hsv = area; kmdl = 1;
    xr[0] = xc; yr[0] = yc; zr[0] = za; zr[1] = za+w; wr = w;
}

void mhrect(double am, double cx, double rx, double cy, double ry,
            double zc, double h)
{   amr = am * 100.;  hsv = h; kmdl = 2; if ((rx==0.) || (ry==0.)) kmdl = -1;
    xr[0] = cx - rx/2.; xr[1] = cx + rx/2.;
    yr[0] = cy - ry/2.; yr[1] = cy + ry/2.; zr[0] = zc;
    dcz = sqrt(abs(rx*ry)) / 100000.;
}

void mprism(double am, double cx, double rx, double cy, double ry,
            double za, double w)
{   amr = am * 100.; kmdl = 3; if ((rx==0.) || (ry==0.)) kmdl = -1;
    xr[0] = cx - rx/2.; xr[1] = cx + rx/2.;
    yr[0] = cy - ry/2.; yr[1] = cy + ry/2.; zr[0] = za; zr[1] = za+w; wr = w;
    dcx = rx / 100000.; dcy = ry / 100000.;
    if (w != 0.) dcz = w / 100000.; else dcz = sqrt(abs(dcx*dcy));
}

double calma(double xp, double yp, double zp)
{   int k, km;
    double x, y, z, xx, yy, zz, rr, r, r3, r5, ss, s4;
    double t1, t2, t3, t4, t5, t6, tt;
    switch (kmdl) {
    case 0:
        t4 = cx-cy; t5 = cy-cz; t6 = cz-cx;
        t1 = t4-t6; t2 = t5-t4; t3 = t6-t5;
        x = xr[0]-xp; xx = x*x;
        y = yr[0]-yp; yy = y*y;
        z = zr[0]-zp; zz = z*z;
        rr = xx+yy+zz; r = sqrt(rr); r5 = rr*rr*r;
        if (r5 == 0.) return(0.);
        tt =  t1*xx + t2*yy + t3*zz + (cxy*x*y + cyz*y*z + czx*z*x)*6.;
        return(amr*hsv*tt/r5);
    case 1:
        t4 = cx-cy; t5 = cy-cz; t6 = cz-cx;
        t1 = t4-t6; t2 = t5-t4; t3 = t6-t5;
        x = xr[0]-xp; xx = x*x; y = yr[0]-yp; yy = y*y;
        ss = xx+yy; z = zr[0]-zp; zz = z*z;
        if (ss == 0.) {
            if (zz == 0.) return(0.);
            if (wr == 0.) tt = (t3/2.) / zz;
            else {
                t1 = zz; z = zr[1]-zp; t2 = z*z;
                if (t2 == 0.) return(0.);
                tt = (t3/2.) * (1./t1 - 1./t2);
            }
        }
        else {
            s4 = ss*ss; rr = ss+zz; r = sqrt(rr); r3 = rr*r;
            t4 = t1*xx + t2*yy + cxy*x*y*6.;
            t5 = t5*yy - t6*xx + cxy*x*y*2.;
            t6 = czx*x + cyz*y;
            t1 = z/r; t2 = zz*z/r3; t3 = 2./r3;
            if (wr == 0.) { t1 = 1.-t1; t2 = 1.-t2; t3 = -t3; }
            else {
                z = zr[1]-zp; zz = z*z; rr = ss+zz; r = sqrt(rr); r3 = rr*r;
                t1 = z/r - t1; t2 = zz*z/r3 - t2; t3 = 2./r3 - t3;
            }
            tt = (t1*t4-t2*t5)/s4 - t3*t6;
        }
        return(amr*hsv*tt);
    case 2:
        t1=0.; t2=0.; t3=0.; t4=0.; t5=0.;
        z = zr[0]-zp; if (z == 0.) z = dcz; zz = z*z;
        for (k=0; k<4; k++) {
            if (k == 2) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; t5 = -t5; }
            x = xr[k%2]-xp; xx = x*x;
            y = yr[k/2]-yp; yy = y*y;
            rr = xx+yy+zz; r = sqrt(rr);
            r3 = r*(xx+zz); t1 = x*y/r3 - t1; t4 = y*z/r3 - t4;
            r3 = r*(yy+zz); t2 = x*y/r3 - t2; t3 = x*z/r3 - t3;
            t5 = 1./r - t5;
        }
        tt = (cz-cx)*t1 + (cz-cy)*t2 - (cyz*t3 + czx*t4 - cxy*t5)*2.;
        return(amr*hsv*tt);
    case 3:
        t1=1.; t2=1.; t3=1.; t4=0.; t5=0.; t6=0.;
        k = 0; km = ((wr==0.)? 4: 8);
        if ((xr[0]==xp) || (xr[1]==xp)) xp -= dcx;
        if ((yr[0]==yp) || (yr[1]==yp)) yp -= dcy;
        if ((zr[0]==zp) || (zr[1]==zp)) zp -= dcz;
        while (k < km) {
            if ((k%4) == 2) {
                t1 = 1./t1; t2 = 1./t2; t3 = 1./t3;
                t4 = -t4; t5 = -t5; t6 = -t6;
            }
            x = xr[k%2]-xp;   xx = x*x;
            y = yr[k/2%2]-yp; yy = y*y;
            z = zr[k/4]-zp;   zz = z*z;
            rr = xx+yy+zz; r = sqrt(rr);
            t1 = ((r-x)/(r+x))/t1; t2 = ((r-y)/(r+y))/t2;
            t3 = ((r-z)/(r+z))/t3;
            t4 = atan((y*z)/(r*x))-t4;
            t5 = atan((z*x)/(r*y))-t5;
            t6 = atan((x*y)/(r*z))-t6;
            k++;
        }
        tt = cyz*log(t1) + czx*log(t2) + cxy*log(t3)
              + cx*t4 + cy*t5 + cz*t6;
        if (km == 8) tt = -tt;
        return(amr*tt);
    default:
        return(0.);
    }
}

/*******************************************/
/*   FORTRAN Language Interface Routines   */
/*******************************************/

#define W   double

void magafd_(float *the, float *phi, float *dip, float *dec)
{   magafd((W)*the, (W)*phi, (W)*dip, (W)*dec); }

void mpoint_(float *am, float *xc, float *yc, float *zc, float *vol)
{   mpoint((W)*am, (W)*xc, (W)*yc, (W)*zc, (W)*vol); }

void mvline_(float *am, float *xc, float *yc, float *za, float *w, float *area)
{   mvline((W)*am, (W)*xc, (W)*yc, (W)*za, (W)*w, (W)*area); }

void mhrect_(float *am, float *cx, float *rx, float *cy, float *ry,
             float *zc, float *h)
{   mhrect((W)*am, (W)*cx, (W)*rx, (W)*cy, (W)*ry, (W)*zc, (W)*h); }

void mprism_(float *am, float *cx, float *rx, float *cy, float *ry,
             float *za, float *w)
{   mprism((W)*am, (W)*cx, (W)*rx, (W)*cy, (W)*ry, (W)*za, (W)*w); }

float calma_(float *xp, float *yp, float *zp)
{   return((float)calma((W)*xp,(W)*yp,(W)*zp)); }