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

#define RAD     (180./3.14159)

static int kmdl;
static double amr, cx, cy, cz, cxy, cyz, czx;
static double xr[2], yr[2], zr[2], wr, hsv;

void magafd(double am, double the, double phi, double dip, double dec)
{   double cthe, cdip, px, py, pz, ex, ey, ez;
    amr = am * 100.;
    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 xc, double yc, double zc, double vol)
{   xr[0] = xc; yr[0] = yc; zr[0] = zc; hsv = vol; kmdl = 0; }

void mvline(double xc, double yc, double zs, double w, double area)
{   xr[0] = xc; yr[0] = yc; zr[0] = zs; zr[1] = zs+w; wr = w;
    hsv = area; kmdl = 1; }

void mhrect(double xs, double xt, double ys, double yt, double zc, double h)
{   xr[0] = xs; xr[1] = xt; yr[0] = ys; yr[1] = yt; zr[0] = zc;
    hsv = h; kmdl = 2; }

void mprism(double xs, double xt, double ys, double yt, double zs, double w)
{   xr[0] = xs; xr[1] = xt; yr[0] = ys; yr[1] = yt; zr[0] = zs; zr[1] = zs+w;
    wr = w; kmdl = 3; }

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;
        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 (wr == 0.) tt = (t3/2.) / zz;
            else {
                t1 = zz; z = zr[1]-zp; t2 = z*z;
                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; 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);
        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)/t3;
            t4 = atan2(x*y,rr+r*z-yy)-t4;
            t5 = atan2(x*y,rr+r*z-xx)-t5;
            t6 = atan((x*y)/(r*z))-t6;
            k++;
        }
        tt =  cyz*log(t1) + czx*log(t2) - cxy*log(t3)*2.
            - cx*t4 - cy*t5 + cz*t6;
        if (km == 8) tt = -tt;
        return(amr*tt);
    }
}

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

#define W   double

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

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

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

void mhrect_(float *xs, float *xt, float *ys, float *yt, float *zc, float *h)
{   mprism((W)*xs, (W)*xt, (W)*ys, (W)*yt, (W)*zc, (W)*h); }

void mprism_(float *xs, float *xt, float *ys, float *yt, float *zs, float *w)
{   mprism((W)*xs, (W)*xt, (W)*ys, (W)*yt, (W)*zs, (W)*w); }

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