Locating the intersection

/* geosi - intersection of two spheres */

geosi(s) double s[][7]; {
int    i, j, k, n, m, m0;
double b1, b2, c1, c2, g, th, dh, ri;
double xx, yy, zz, rr;
double w[3], u[3], o[3][3];
n=4096;
n=2048;
g=128.0;
g=64.0;
g=53.0;
dh=3.14159/((double)n);
for (i=0; i<2; i++) 
for (j=0; j<2; j++) {
  if (i==j) continue;
  th=0.0;
  ri=s[i][3];
/*  twopv(&s[i][0]); */
/*  twopv(&s[i][4]); */
  spheu(o,&s[i][4]);
/*  twopm(o); */
  sphrv(w,ri,th,g*th);
  sphap(u,o,w,&s[i][0]);
  c1=-1.0;
  c2=-1.0;
  m0=3;
  pltms(u[0],u[2],0);
  for (k=0; k<n; k++) {
    th+=dh;
    sphrv(w,ri,th,g*th);
    sphap(w,o,w,&s[i][0]);
    xx=w[0]-s[j][0]; xx*=xx;
    yy=w[1]-s[j][1]; yy*=yy;
    zz=w[2]-s[j][2]; zz*=zz;
    rr=s[j][3]; rr*=rr;
    b2=xx+zz-rr;
    b1=b2+yy;
    m=3;
    if      (s[i][1]<w[1]) m0=4;
    else if (b1<0.0) {m=1; m0=1;}
    else if (w[1]<=s[j][1]) m0=3;
    else if (b2>=0.0) m0=3;
    else    m=2;
    switch (m0) {
      case 1: pltil(u[0],u[2],c1,w[0],w[2],b1,i+1); break;
      case 2: pltil(u[0],u[2],c2,w[0],w[2],b2,i+1); break;
      case 3: pltil(u[0],u[2],1.0,w[0],w[2],1.0,i+1); break;
      case 4: pltil(u[0],u[2],-1.0,w[0],w[2],-1.0,i+1); break;
      case 5: pltms(w[0],w[2],i+1); break;
      default: break;
      }

    u[0]=w[0]; u[1]=w[1]; u[2]=w[2];
    c1=b1;
    c2=b2;
    m0=m;
    }  /* end for k */
  }    /* end for j */
}