/* Subroutine to interpolate the position of a contour which is known */
/* to be next to ix,iy in the direction ixg,iyg. The unscaled distance */
/* along ixg,iyg is returned as dist */

#include "plplot.h"

void plccal(pts,nx,ny,zlev,ix,iy,ixg,iyg,dist)
int nx,ny,ix,iy,ixg,iyg;
double zlev,*dist;
FLOAT **pts;
{
      int ia,ib;
      double dbot, dtop, pmid, qmid, zmid;

      ia = ix+ixg;
      ib = iy+iyg;
      if (ixg == 0 || iyg == 0) {
        dtop = zlev - pts[ix-1][iy-1];
        dbot = pts[ia-1][ib-1] - pts[ix-1][iy-1];
        *dist = 0.0;
        if (dbot != 0.0) *dist = dtop/dbot;
      }
      else {
        pmid = pts[ix-1][iy-1] + pts[ia-1][ib-1];
        qmid = pts[ix-1][ib-1] + pts[ia-1][iy-1];
        zmid = (pmid+qmid)/4.0;
        if (zmid >= zlev) {
          dtop = zlev - pts[ix-1][iy-1];
          dbot = zmid - pts[ix-1][iy-1];
          *dist = 0.0;
          if (dbot != 0.0) *dist = 0.5*dtop/dbot;
        }
        else {
          dtop = zlev-zmid;
          dbot = pts[ia-1][ib-1] - zmid;
          *dist = 0.5;
          if (dbot != 0.0) *dist = 0.5 + 0.5*dtop/dbot;
        }
      }
}