/* 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;
float *pts,zlev,*dist;
{
      int ia,ib;
      float dbot, dtop, pmid, qmid, zmid;

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