#include <stdlib.h> 
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include "xpplim.h"
#include "struct.h"

#define MAX_LEN_SBOX 25
#define DING ping
#define MAX_NULL 10000

extern int SuppressBounds;
extern GRAPH *MyGraph;

int NullStyle=0; /* 1 is with little vertical/horizontal lines */
extern (*rhs)();

double atof();
extern int DRight,DLeft,DTop,DBottom;
extern int STORFLAG;
extern double last_ic[MAXODE];

extern double DELTA_T,TEND,TRANS;
extern int PaperWhite,DCURY;
int XNullColor=2,YNullColor=7;
int NULL_HERE,num_x_n,num_y_n,num_index,
	null_ix,null_iy,WHICH_CRV;
float null_dist,*X_n,*Y_n,*saver,*NTop,*NBot;
extern int NMESH,NODE,NJMP,NMarkov,FIX_VAR,NEQ;
float fnull();
int DF_GRID=10,DF_FLAG=0,DF_IX=-1,DF_IY=-1;
int DFIELD_TYPE=0;

typedef struct {
		float x,y,z;
		} Pt;

typedef struct nclines {
                float *xn,*yn;
                int nmx,nmy;
                int n_ix,n_iy;
                struct nclines *n,*p;
}  NCLINES;


RANGE_INFO ncrange;  

NCLINES *ncperm;
int n_nstore=0;
int ncline_cnt;

froz_cline_stuff_com(int i)
{
  int delay=200;
  if(n_nstore==0)start_ncline();
  switch(i){
  case 0:
    if(NULL_HERE==0)return;
    add_froz_cline(X_n,num_x_n,null_ix,Y_n,num_y_n,null_iy);
    break;
  case 1:
    clear_froz_cline();
    break;
  case 3:
    new_int("Delay (msec)",&delay);
    if(delay<=0)delay=0;
    redraw_froz_cline(delay);
    break;
  case 2:
    do_range_clines();
    break;
  }
 }

	
		 

do_range_clines()
{
  static char *n[]={"*2Range parameter","Steps","Low","High"};
  char values[4][MAX_LEN_SBOX];
  int status,i;
  double z,dz,zold;
  float xmin,xmax,y_tp,y_bot;
  int col1=XNullColor,col2=YNullColor;
  int course=NMESH;
  /* if(PaperWhite){
    col1=1;
    col2=9;
    } */
  sprintf(values[0],"%s",ncrange.rv);
  sprintf(values[1],"%d",ncrange.nstep);
  sprintf(values[2],"%g",ncrange.xlo);
  sprintf(values[3],"%g",ncrange.xhi);
  status=do_string_box(4,4,1,"Range Clines",n,values,45);
  if(status!=0){
    strcpy(ncrange.rv,values[0]);
    ncrange.nstep=atoi(values[1]);
    ncrange.xlo=atof(values[2]);
    ncrange.xhi=atof(values[3]);
    if(ncrange.nstep<=0)return;
    dz=(ncrange.xhi-ncrange.xlo)/(double)ncrange.nstep;
    if(dz<=0.0)return;
    get_val(ncrange.rv,&zold);
    
    for(i=NODE;i<NODE+NMarkov;i++)set_ivar(i+1+FIX_VAR,last_ic[i]);
    xmin=(float)MyGraph->xmin;
    xmax=(float)MyGraph->xmax;
    y_tp=(float)MyGraph->ymax;
    y_bot=(float)MyGraph->ymin;
    null_ix=MyGraph->xv[0];
    null_iy=MyGraph->yv[0];
    
    
    for(i=0;i<=ncrange.nstep;i++){
      z=(double)i*dz+ncrange.xlo;
      set_val(ncrange.rv,z);
      if(NULL_HERE==0)
	{
	  if((X_n=(float *)malloc(4*MAX_NULL*sizeof(float)))!=NULL
	     && (Y_n=(float *)malloc(4*MAX_NULL*sizeof(float)))!=NULL)
	    
	    
	    NULL_HERE=1;
	  NTop=(float *)malloc((course+1)*sizeof(float));
	  NBot=(float *)malloc((course+1)*sizeof(float));
	  if(NTop==NULL||NBot==NULL)NULL_HERE=0;
	}
      else {
	free(NTop);
	free(NBot);
	NTop=(float *)malloc((course+1)*sizeof(float));
	NBot=(float *)malloc((course+1)*sizeof(float));
	if(NTop==NULL||NBot==NULL){NULL_HERE=0;
	return;}
      }
      
      WHICH_CRV=null_ix;
      set_linestyle(col1);
      new_nullcline(course,xmin,y_bot,xmax,y_tp,X_n,&num_x_n);
      
      
      WHICH_CRV=null_iy;
      set_linestyle(col2);
      new_nullcline(course,xmin,y_bot,xmax,y_tp,Y_n,&num_y_n);
      add_froz_cline(X_n,num_x_n,null_ix,Y_n,num_y_n,null_iy);
    }
    set_val(ncrange.rv,zold);
  }
  
}

start_ncline()
{
  n_nstore=1;
  ncperm=(NCLINES *)malloc(sizeof(NCLINES));
  ncperm->p=NULL;
  ncperm->n=NULL;
  ncperm->nmx=0;
  ncperm->nmy=0;
  ncperm->n_ix=-5;
  ncperm->n_iy=-5;
  ncrange.xlo=0;
  ncrange.xhi=1;
  ncrange.nstep=10;
  sprintf(ncrange.rv," ");
}

clear_froz_cline()
{
  NCLINES *z,*znew;
  z=ncperm;
  while(z->n!=NULL)
    z=z->n;
  /*  this is the bottom but there is nothing here that has been stored   */
  
  znew=z->p;
  if(znew==NULL)return;
  free(z);
  z=znew;
  /* now we are deleting everything */
   while(z->p !=NULL){
      znew=z->p;
      z->n=NULL;
      z->p=NULL;
      free(z->xn);
      free(z->yn);
      free(z);
      z=znew;
  }
  if(ncperm->nmx>0){
    free(ncperm->xn);
    ncperm->nmx=0;
  }
  if(ncperm->nmy>0){
    free(ncperm->yn);
    ncperm->nmy=0;
  }
  ncperm->n=NULL;
  n_nstore=1;
  ncline_cnt=0;
}

get_nullcline_floats(float **v,int *n,int who,int type) /* type=0,1 */
{
  NCLINES *z;
  int i;
  if(who<0){
    if(type==0){
      *v=X_n;
      *n=num_x_n;
    }
    else {
      *v=Y_n;
      *n=num_y_n;
    }
    if(v==NULL)return 1;
    return 0;
  }
  if(who>ncline_cnt||n_nstore==0)return 1;
   z=ncperm;
   for(i=0;i<who;i++)
     z=z->n;
   if(z==NULL)return 1;
   if(type==0){
      *v=z->xn;
      *n=z->nmx;
    }
    else {
      *v=z->yn;
      *n=z->nmy;
    }
    if(v==NULL)return 1;
    return 0;
}   

redraw_froz_cline(flag)
     int flag;
{
  NCLINES *z;
  int col1=XNullColor,col2=YNullColor;
  /* if(PaperWhite){
    col1=1;
    col2=9;
    } */
  if(n_nstore==0)return;
  z=ncperm;
  while(1){
    if(z==NULL||(z->nmx==0&&z->nmy==0))return;
    
  /*  printf(" %d %d  %d %d  %d \n",
	   MyGraph->xv[0],z->n_ix, &MyGraph->yv[0],z->n_iy ,
	  MyGraph->ThreeDFlag==0); */
    if(MyGraph->xv[0]==z->n_ix&&MyGraph->yv[0]==z->n_iy
       &&MyGraph->ThreeDFlag==0)
      {
	if(flag>0){
	  waitasec(flag);
	  clr_scrn();
	}
	set_linestyle(col1);
	restor_null(z->xn,z->nmx,1);
	set_linestyle(col2);
	restor_null(z->yn,z->nmy,2);
	if(flag>0)
	  FlushDisplay();
      }
    	z=z->n;
	if(z==NULL)break;

    
  }
}

add_froz_cline(xn,nmx,n_ix,yn,nmy,n_iy)
     float *xn,*yn;
     int nmx,nmy,n_ix,n_iy;
{
  NCLINES *z,*znew;
  int i;
  z=ncperm;
  /* move to end */
  while(z->n!=NULL){
    z=(z->n); 
  }
  z->xn=(float *)malloc(4*nmx*sizeof(float));
  for(i=0;i<4*nmx;i++)
    z->xn[i]=xn[i];
  z->yn=(float *)malloc(4*nmy*sizeof(float));
  for(i=0;i<4*nmy;i++)
    z->yn[i]=yn[i]; 
  z->nmx=nmx;
  z->nmy=nmy;
  z->n_ix=n_ix;
  z->n_iy=n_iy;
  z->n=(NCLINES *)malloc(sizeof(NCLINES));
  znew=z->n;
  znew->n=NULL;
  znew->p=z;
  znew->nmx=0;
  znew->nmy=0;
  znew->n_ix=-5;
  znew->n_iy=-5;
  ncline_cnt++;
}

                
get_max_dfield(y,ydot,u0,v0,du,dv,n,inx,iny,mdf)
     double *y,*ydot,du,dv,u0,v0,*mdf;
     int n,inx,iny;
{
  int i,j;
  double amp,dxp,dyp;
  *mdf=0.0;
  for(i=0;i<=n;i++){
    y[inx]=u0+du*i;
    for(j=0;j<=n;j++){
      y[iny]=v0+dv*j;
      rhs(0.0,y,ydot,NODE);
      extra(y,0.0,NODE,NEQ);
      scale_dxdy(ydot[inx],ydot[iny],&dxp,&dyp);
      amp=hypot(dxp,dyp);
      if(amp>*mdf)*mdf=amp;
    }
  }
}
/*  all the nifty 2D stuff here    */

redraw_dfield()
{
  int i,j,start,k;
  int inx=MyGraph->xv[0]-1;
  int iny=MyGraph->yv[0]-1;
  double y[MAXODE],ydot[MAXODE],xv1,xv2;
  float v1[MAXODE],v2[MAXODE];
  double scx=MyGraph->yhi-MyGraph->ylo;
  double scy=MyGraph->xhi-MyGraph->xlo;
  double amp,mdf;
  double t;
  double du,dv,u0,v0,dxp,dyp,dz,dup,dvp;
  double oldtrans=TRANS;

  int grid=DF_GRID;
  
  
  if(DF_FLAG==0||
     MyGraph->TimeFlag||MyGraph->xv[0]==MyGraph->yv[0]||MyGraph->ThreeDFlag
     || DF_IX!=MyGraph->xv[0]||DF_IY!=MyGraph->yv[0])
    return;
  
  du=(MyGraph->xhi-MyGraph->xlo)/(double)grid;
  dv=(MyGraph->yhi-MyGraph->ylo)/(double)grid;
  
  dup =(double)(DRight-DLeft)/(double)grid;
  dvp=(double)(DTop-DBottom)/(double)grid; 
  dz=hypot(dup,dvp)*(.25+.75*DFIELD_TYPE);
  u0=MyGraph->xlo;
  v0=MyGraph->ylo;
  set_linestyle(MyGraph->color[0]);
  get_ic(2,y);
  get_max_dfield(y,ydot,u0,v0,du,dv,grid,inx,iny,&mdf);
  for(i=0;i<=grid;i++){
    y[inx]=u0+du*i;
    for(j=0;j<=grid;j++){
      y[iny]=v0+dv*j;
      rhs(0.0,y,ydot,NODE);
      extra(y,0.0,NODE,NEQ);
      if(MyGraph->ColorFlag||DF_FLAG==2){
	v1[0]=0.0;
	v2[0]=0.0;
	for(k=0;k<NEQ;k++){
	  v1[k+1]=(float)y[k];
	  v2[k+1]=v1[k+1]+(float)ydot[k];
	}
	comp_color(v1,v2,NODE,1.0);
      }
      if(DF_FLAG==1){
	scale_dxdy(ydot[inx],ydot[iny],&dxp,&dyp);
	if(DFIELD_TYPE==1)
	  {
	    ydot[inx]/=mdf;
	    ydot[iny]/=mdf;
	    }
	else{
	  amp=hypot(dxp,dyp);
	  if(amp!=0.0){
	    ydot[inx]/=amp;
	    ydot[iny]/=amp;
	  }
	}
	xv1=y[inx]+ydot[inx]*dz;
	xv2=y[iny]+ydot[iny]*dz;
        bead_abs((float)xv1,(float)xv2);
	line_abs((float)y[inx],(float)y[iny],(float)xv1,(float)xv2);
      }
      if(DF_FLAG==2&&j>0&&i<grid){
	frect_abs((float)y[inx],(float)y[iny],(float)du,(float)dv);
      }
    }
  }
}

direct_field_com(int c)
{
  
  int i,j,start,k;
  int inx=MyGraph->xv[0]-1;
  int iny=MyGraph->yv[0]-1;
  double y[MAXODE],ydot[MAXODE],xv1,xv2;
  double dtold=DELTA_T;
  float v1[MAXODE],v2[MAXODE];
  double scx=MyGraph->yhi-MyGraph->ylo;
  double scy=MyGraph->xhi-MyGraph->xlo;
  double amp,mdf;
  double t;
  double du,dv,u0,v0,dxp,dyp,dz,dup,dvp;
  double oldtrans=TRANS;
  
  
  int grid=DF_GRID;
  
  
  if(MyGraph->TimeFlag||MyGraph->xv[0]==MyGraph->yv[0]||MyGraph->ThreeDFlag)
    return;

  if(c==2){
    DF_FLAG=0;
    return;
  }
  if(c==0)DFIELD_TYPE=1;
  if(c==4)DFIELD_TYPE=0;
  new_int("Grid:",&grid);
  if(grid<=1)return;
  
  du=(MyGraph->xhi-MyGraph->xlo)/(double)grid;
  dv=(MyGraph->yhi-MyGraph->ylo)/(double)grid;
  
  dup =(double)(DRight-DLeft)/(double)grid;
  dvp=(double)(DTop-DBottom)/(double)grid; 
  dz=hypot(dup,dvp)*(.25+.75*DFIELD_TYPE) ;
  u0=MyGraph->xlo;
  v0=MyGraph->ylo;
  set_linestyle(MyGraph->color[0]);
  if(c!=1){
    DF_GRID=grid;
    DF_FLAG=1;
    if(c==3)
      DF_FLAG=2;
    DF_IX=inx+1;
    DF_IY=iny+1;
    get_ic(2,y);
     get_max_dfield(y,ydot,u0,v0,du,dv,grid,inx,iny,&mdf);
   for(i=0;i<=grid;i++){
     y[inx]=u0+du*i;
     for(j=0;j<=grid;j++){
       y[iny]=v0+dv*j;
       rhs(0.0,y,ydot,NODE);
       extra(y,0.0,NODE,NEQ);
       if(MyGraph->ColorFlag||DF_FLAG==2){
	 v1[0]=0.0;
         v2[0]=0.0;
	 for(k=0;k<NEQ;k++){
	   v1[k+1]=(float)y[k];
	   v2[k+1]=v1[k+1]+(float)ydot[k];
	 }
	 comp_color(v1,v2,NODE,1.0);
       }
       if(DF_FLAG==1){
	 scale_dxdy(ydot[inx],ydot[iny],&dxp,&dyp);
	 if(DFIELD_TYPE==0){
	   amp=hypot(dxp,dyp);
	   if(amp!=0.0){
	     ydot[inx]/=amp;
	     ydot[iny]/=amp;
	   }
	 }
	 else {
	   ydot[inx]/=mdf;
	   ydot[iny]/=mdf;
	 }
	   xv1=y[inx]+ydot[inx]*dz;
	 xv2=y[iny]+ydot[iny]*dz;
         bead_abs((float)xv1,(float)xv2);
	 line_abs((float)y[inx],(float)y[iny],(float)xv1,(float)xv2);
       }
       if(DF_FLAG==2&&j>0&&i<grid){
	 frect_abs((float)y[inx],(float)y[iny],(float)du,(float)dv);
       }
       
     }
   }
   TRANS=oldtrans;
   return;
  }
  STORFLAG=0;
   SuppressBounds=1;
   for(k=0;k<2;k++){
     for(i=0;i<=grid;i++)
       for(j=0;j<=grid;j++)
	 {
	 get_ic(2,y);
	 y[inx]=u0+du*i;
	 y[iny]=v0+dv*j;
	 t=0.0;
	 start=1;
	 if(integrate(&t,y,TEND,DELTA_T,1,NJMP,&start)==1){
	   TRANS=oldtrans;
	   return;
	   STORFLAG=1;
	 }
	 }
     DELTA_T=-DELTA_T;
   }
   SuppressBounds=0;
   DELTA_T=dtold;
}

/* animated nullclines stuff   
   added Aug 31 97 
   just redraws them
   It will allow you to either freeze a range of them
   or just one at a time

   freeze   -   store the current set
   range_freeze - compute over some range of parameters
   clear - delete all but the current set
   animate - replay all frozen ones (not current set )
   */ 	 
       
       
     
     
save_the_nullclines()
{
  FILE *fp;
  char filename[256];
  if(NULL_HERE==0)return;
  sprintf(filename,"nc.dat");
  ping();
  if(!file_selector("Save nullclines",filename,"*.dat"))return;
  fp=fopen(filename,"w");
  if(fp==NULL){
    err_msg("Cant open file!");
    return;
  }
  dump_clines(fp,X_n,num_x_n,Y_n,num_y_n);
  fclose(fp);
}

restore_nullclines()
{
 int col1=XNullColor,col2=YNullColor;
 /* if(PaperWhite){
   col1=1;
   col2=9;
   } */
 if(NULL_HERE==0)return;
 if(MyGraph->xv[0]==null_ix&&MyGraph->yv[0]==null_iy&&MyGraph->ThreeDFlag==0)
   {
    set_linestyle(col1);
    restor_null(X_n,num_x_n,1);
    set_linestyle(col2);
    restor_null(Y_n,num_y_n,2);
  }
 redraw_froz_cline(0);
}

dump_clines(fp,x,nx,y,ny) /* gnuplot format */
     FILE *fp;
     float *x,*y;
     int nx,ny;
{
    int i;
    fprintf(fp,"# X-nullcline\n");
    for(i=0;i<nx-1;i++){
      fprintf(fp,"%g %g 1 \n",x[4*i],x[4*i+1]);
      fprintf(fp,"%g %g 1 \n",x[4*i+2],x[4*i+3]);
      fprintf(fp,"\n");
    }
    fprintf(fp,"\n# Y-nullcline\n");
     for(i=0;i<ny-1;i++){
      fprintf(fp,"%g %g 2 \n",y[4*i],y[4*i+1]);
      fprintf(fp,"%g %g 2 \n",y[4*i+2],y[4*i+3]);
      fprintf(fp,"\n");
    }


}
dump_clines_old(fp,x,nx,y,ny)
     FILE *fp;
     float *x,*y;
     int nx,ny;
{
    int i,ix,iy;
    int n;
    n=nx;
    if(n<ny)n=ny;
    for(i=0;i<n;i++){
      if(i>=nx)
	ix=nx-1;
      else 
	ix=i;
      if(i>=ny)
	iy=ny-1;
      else
	iy=i;
      fprintf(fp,"%g %g %g %g \n",x[4*ix],x[4*ix+1],y[4*iy],y[4*iy+1]);
      fprintf(fp,"%g %g %g %g \n",x[4*ix+2],x[4*ix+3],y[4*iy+2],y[4*iy+3]);
      
      
    }

}
restor_null(v,n,d) /* d=1 for x and 2 for y  */
     float *v;
     int n,d;
{
  
  int i,i4;
  float xm,ym;
  int x1,y1;
  for(i=0;i<n;i++)
    {
      i4=4*i;
      line_abs(v[i4],v[i4+1],v[i4+2],v[i4+3]);
      if(NullStyle==1){
	xm=.5*(v[i4]+v[i4+2]);
	ym=.5*(v[i4+1]+v[i4+3]);
	scale_to_screen(xm,ym,&x1,&y1);
	switch(d){
	case 1: 
	  line(x1,y1-4,x1,y1+4);
	  
	  break;
	case 2:
	  line(x1-4,y1,x1+4,y1);
	  break;
	    }
      }
    }
}


new_clines_com(int c)
{
  int course=NMESH,i;
  float xmin,xmax,y_tp,y_bot;
  int col1=XNullColor,col2=YNullColor;
  
  if(MyGraph->ThreeDFlag||MyGraph->TimeFlag||MyGraph->xv[0]==MyGraph->yv[0])return;

  if(c==1){
    restore_nullclines();
    return;
  }
  if(c==2){
    MyGraph->Nullrestore=1;
   return;
  }
  if(c==3){
    MyGraph->Nullrestore=0;
    return;
  }
  if(c==4){
    froz_cline_stuff();
    return;
  }
  if(c==5){
    save_the_nullclines();
    return;
  }
  if(c==0){
    for(i=NODE;i<NODE+NMarkov;i++)set_ivar(i+1+FIX_VAR,last_ic[i]);
    xmin=(float)MyGraph->xmin;
    xmax=(float)MyGraph->xmax;
    y_tp=(float)MyGraph->ymax;
    y_bot=(float)MyGraph->ymin;
  null_ix=MyGraph->xv[0];
  null_iy=MyGraph->yv[0];
  if(NULL_HERE==0)
    {
      if((X_n=(float *)malloc(4*MAX_NULL*sizeof(float)))!=NULL
	 && (Y_n=(float *)malloc(4*MAX_NULL*sizeof(float)))!=NULL)
	
	
	NULL_HERE=1;
      NTop=(float *)malloc((course+1)*sizeof(float));
      NBot=(float *)malloc((course+1)*sizeof(float));
      if(NTop==NULL||NBot==NULL)NULL_HERE=0;
    }
  else {
    free(NTop);
    free(NBot);
    NTop=(float *)malloc((course+1)*sizeof(float));
   NBot=(float *)malloc((course+1)*sizeof(float));
   if(NTop==NULL||NBot==NULL){NULL_HERE=0;
   return;}
  }
  
  WHICH_CRV=null_ix;
  set_linestyle(col1);
  new_nullcline(course,xmin,y_bot,xmax,y_tp,X_n,&num_x_n);
  ping();
  
  WHICH_CRV=null_iy;
  set_linestyle(col2);
  new_nullcline(course,xmin,y_bot,xmax,y_tp,Y_n,&num_y_n);
  ping();
  
  }
}


new_nullcline(course,xlo,ylo,xhi,yhi,stor,npts)
     int course;
     float xlo,ylo,xhi,yhi;
     int *npts;
     float *stor;
{
 num_index=0;
 saver=stor;
 do_cline(course,xlo,ylo,xhi,yhi);
 *npts=num_index;
}



stor_null(x1,y1,x2,y2)
float x1,y1,x2,y2;
{
 int i;
 if(num_index>=MAX_NULL)return;
 i=4*num_index;
 saver[i]=x1;
 saver[i+1]=y1;
 saver[i+2]=x2;
 saver[i+3]=y2;
 num_index++;
} 

float fnull( x, y)
 float x,y;
 {
  double y1[MAXODE],ydot[MAXODE];
  int i;
  for(i=0;i<NODE;i++)y1[i]=last_ic[i];
 
  y1[null_ix-1]=(double)x;
  y1[null_iy-1]=(double)y;
  rhs(0.0,y1,ydot,NODE);
  /*  printf(" %f  %f %f \n ", x,y,ydot[WHICH_CRV-1]); */
  return((float)ydot[WHICH_CRV-1]);
 }


interpolate(p1,p2,z,x,y)
 Pt p1,p2;
 float z,*x,*y;
{
 float scale;
  if(p1.z==p2.z)return(0);
  scale=(z-p1.z)/(p2.z-p1.z);
  *x=p1.x+scale*(p2.x-p1.x);
  *y=p1.y+scale*(p2.y-p1.y);
   return(1);
 }

quad_contour(p1,p2,p3,p4)
Pt p1,p2,p3,p4;
{
 float x[4],y[4];
 int count=0;
 if(p1.z*p2.z<=0.0)
   if(interpolate(p1,p2,0.0,&x[count],&y[count]))count++;
 if(p2.z*p3.z<=0.0)
   if(interpolate(p3,p2,0.0,&x[count],&y[count]))count++;
 if(p3.z*p4.z<=0.0)
   if(interpolate(p3,p4,0.0,&x[count],&y[count]))count++;
 if(p1.z*p4.z<=0.0)
   if(interpolate(p1,p4,0.0,&x[count],&y[count]))count++;


 if(count==2){
   line_abs(x[0],y[0],x[1],y[1]);
   stor_null(x[0],y[0],x[1],y[1]);
 }
 

}


triangle_contour(p1,p2,p3)

Pt p1,p2,p3;
{
 float x[3],y[3];
 int count=0;
 if(p1.z*p2.z<=0.0)
 /* if(((0.0<=p1.z)&&(0.0>=p2.z))||
	((0.0>=p1.z)&&(0.0<=p2.z))) */
	if(interpolate(p1,p2,0.0,&x[count],&y[count]))count++;
if( p1.z*p3.z<=0.0)
/*  if(((0.0<=p1.z)&&(0.0>=p3.z))||
	((0.0>=p1.z)&&(0.0<=p3.z))) */

	if(interpolate(p1,p3,0.0,&x[count],&y[count]))count++;
if(p2.z*p3.z<=0.0) 
  /* if(((0.0<=p3.z)&&(0.0>=p2.z))||
	((0.0>=p3.z)&&(0.0<=p2.z))) */
	if(interpolate(p3,p2,0.0,&x[count],&y[count]))count++;
 
 if(count==2){
   line_abs(x[0],y[0],x[1],y[1]);
   stor_null(x[0],y[0],x[1],y[1]);
 }
 

 }




do_cline(ngrid,x1,y1,x2,y2)
int ngrid;
float x1,y1,x2,y2;
{
 float dx=(x2-x1)/(float)ngrid;
 float dy=(y2-y1)/(float)ngrid;
 float x,y,z;
 Pt p[5];
 char esc;
 int i,j,cwidth;
 int nx=ngrid+1;
 int ny=ngrid+1;
 cwidth=get_command_width();
 y=y2;
 for(i=0;i<nx;i++){
   x=x1+i*dx;
   NBot[i]=fnull(x,y);
 }

 for(j=1;j<ny;j++){
  plot_command(ny,j,cwidth);
  esc=my_abort();
  if(esc==27)return;
   y=y2-j*dy;
   NTop[0]=NBot[0];
   NBot[0]=fnull(x1,y);
   for(i=1;i<nx;i++){
     x=x1+i*dx;
     NTop[i]=NBot[i];
     NBot[i]=fnull(x,y);
     p[0].x=x-dx;
     p[0].y=y+dy;
     p[0].z=NTop[i-1];
     p[1].x=x;
     p[1].y=y+dy;
     p[1].z=NTop[i];
     p[3].x=x-dx;
     p[3].y=y;
     p[3].z=NBot[i-1];
     p[2].x=x;
     p[2].y=y;
     p[2].z=NBot[i];
 /*      Uncomment for triangle contour   
      p[4].x=.25*(p[0].x+p[1].x+p[2].x+p[3].x);	
     p[4].y=.25*(p[0].y+p[1].y+p[2].y+p[3].y);
     p[4].z=.25*(p[0].z+p[1].z+p[2].z+p[3].z); 

      
     triangle_contour(p[0],p[1],p[4]);
     triangle_contour(p[1],p[4],p[2]);
     triangle_contour(p[4],p[3],p[2]);
     triangle_contour(p[0],p[4],p[3]); */
 /*   Uncomment for quad contour     */
     quad_contour(p[0],p[1],p[2],p[3]); 
     FlushDisplay();
   }
 }

}