ketcindy 20180614.2

git-svn-id: svn://tug.org/texlive/trunk@48014 c570f23f-e606-0410-a88d-b1316a301751

5 files changed, 1216 insertions, 943 deletions

diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
index 27a5434fe41..84de58498f3 100644
--- a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
@@ -108,6 +108,15 @@ void pull6(int num, double ptL[][6],double pt[6]){
   pt[3]=ptL[num][3]; pt[4]=ptL[num][4]; pt[5]=ptL[num][5]; 
 }
 
+int add1(double gdata[], double x){
+  int nall;
+  nall=length1(gdata);
+  nall++;
+  gdata[nall]=x;
+  gdata[0]=nall;
+  return nall;
+}
+
 int add2(double gdata[][2], double x, double y){
   int nall;
   nall=push2(x, y, gdata);
@@ -189,6 +198,28 @@ int appendd3(int level, int from,int upto, double gdata[][3],double mat[][3]){
   return nall+np;
 }
 
+int appendd6(int level, int from,int upto, double gdata[][6],double mat[][6]){
+  int j, np, nall=length6(mat);
+  np=upto-from+1;
+  if(np<0){np=0;}
+  if(np>0){
+    if(level>0){
+      if(nall>0){
+        nall=push6(Inf,level,0,0,0,0,mat);
+      }
+    }
+    for(j=1; j<=np; j++){
+      mat[nall+j][0]=gdata[from+j-1][0];
+      mat[nall+j][1]=gdata[from+j-1][1];
+      mat[nall+j][2]=gdata[from+j-1][2];
+      mat[nall+j][3]=gdata[from+j-1][3];
+      mat[nall+j][4]=gdata[from+j-1][4];
+      mat[nall+j][5]=gdata[from+j-1][5];    }
+  }
+  mat[0][0]=nall+np;
+  return nall+np;
+}
+
 int prependd2(int from,int upto, double gdata[][2],double mat[][2]){
   int j,np, nall=length2(mat);
   if(from<=upto){
@@ -216,7 +247,38 @@ int prependd2(int from,int upto, double gdata[][2],double mat[][2]){
   return nall+np;
 }
 
-int appendd1(double numvec[], double num){
+int prependd3(int from,int upto, double gdata[][3],double mat[][3]){
+  int j,np, nall=length3(mat);
+  if(from<=upto){
+    np=upto-from+1;
+    for(j=nall; j>=1; j--){
+      mat[np+j][0]=mat[j][0];
+      mat[np+j][1]=mat[j][1];
+      mat[np+j][2]=mat[j][2];
+    }
+    for(j=1; j<=np; j++){
+      mat[j][0]=gdata[from+j-1][0];
+      mat[j][1]=gdata[from+j-1][1];
+      mat[j][2]=gdata[from+j-1][2];
+    }
+  }else{
+    np=from-upto+1;
+    for(j=nall; j>=1; j--){
+      mat[np+j][0]=mat[j][0];
+      mat[np+j][1]=mat[j][1];
+      mat[np+j][2]=mat[j][2];
+    }
+    for(j=1; j<=np; j++){
+      mat[nall+j][0]=gdata[from-j+1][0];
+      mat[nall+j][1]=gdata[from-j+1][1];
+      mat[nall+j][2]=gdata[from-j+1][2];
+    }
+  }
+  mat[0][0]=nall+np;
+  return nall+np;
+}
+
+int appendd1(double num,double numvec[]){
   int n;
   n=floor(numvec[0]+0.5);
   numvec[n+1]=num;
@@ -224,7 +286,7 @@ int appendd1(double numvec[], double num){
   return n+1;
 }
 
-int appendi1(int numvec[],int num){
+int appendi1(int num,int numvec[]){
   int n;
   n=numvec[0];
   numvec[n+1]=num;
@@ -232,6 +294,17 @@ int appendi1(int numvec[],int num){
   return n+1;
 }
 
+int prependd1(double num, double numvec[]){
+  int i,n;
+  n=numvec[0];
+  for(i=n; i>=1; i--){
+    numvec[i+1]=numvec[i];
+  }
+  numvec[1]=num;
+  numvec[0]=n+1;
+  return n+1;
+}
+
 int prependi1(int num, int numvec[]){
   int i,n;
   n=numvec[0];
@@ -259,6 +332,19 @@ int removei2(int nrmv,int mat[][2]){
   }
 }
 
+void removed3(int nrmv,double mat[][3]){
+  int j,nall=length3(mat);
+  if(nall>0){
+    for(j=nrmv; j<=nall-1; j++){
+      mat[j][0]=mat[j+1][0];
+      mat[j][1]=mat[j+1][1];
+      mat[j][2]=mat[j+1][2];
+    }
+    mat[0][0]=nall-1;
+    return;
+  }
+}
+
 int dataindexd1(double out[],int din[][2]){
   int n1=1, j, ctr=0, nall=floor(out[0]+0.5);
   if(nall==0){
@@ -349,6 +435,29 @@ int dataindexd4(int level, double out[][4],int din[][2]){
   return ctr;
 }
 
+int dataindexd6(int level, double out[][6],int din[][2]){
+  int n1=1, j, ctr=0, nall=length6(out);
+  if(nall==0){
+    din[0][0]=0;
+    return 0;
+  }
+  for(j=1; j<=nall; j++){
+    if(out[j][0]>Inf-1){
+      if((level==0) || (out[j][1]==level)){
+        ctr++;
+        din[ctr][0]=n1;
+        din[ctr][1]=j-1;
+        n1=j+1;
+      }
+    }
+  }
+  ctr++;
+  din[ctr][0]=n1;
+  din[ctr][1]=nall;
+  din[0][0]=ctr;
+  return ctr;
+}
+
 int dataindexi1(int out[],int din[][2]){
   int n1=1, j, ctr=0, nall=out[0], ninfp=0;
   if(nall==0){
@@ -379,6 +488,102 @@ int dataindexi1(int out[],int din[][2]){
   return ctr;
 }
 
+void replacelistd1(int ii,double out[],double rep[]){
+  int din[DsizeS][2],k,j,js,je, nall=length1(out),nr=length1(rep);
+  double tmpd1[DsizeL];
+  dataindexd1(out,din);
+  tmpd1[0]=0;
+  for(k=1;k<=din[0][0];k++){
+    if(k!=1){
+      appendd1(Inf,tmpd1);
+    }
+    if(k!=ii){
+      js=din[k][0]; je=din[k][1];
+      for(j=js;j<=je;j++){
+        appendd1(out[j],tmpd1);
+      }
+    }else{
+      for(j=1;j<=nr;j++){
+        appendd1(rep[j],tmpd1);
+      }
+    }
+  }
+  nall=length1(tmpd1);
+  out[0]=0;
+  for(j=1;j<=nall;j++){
+    appendd1(tmpd1[j],out);
+  }
+}
+
+void replacelistd2(int level,int ii,double out[][2],double rep[][2]){//180419
+  int din[DsizeS][2],k,js,je, nall=length2(out),nr=length2(rep);
+  double tmpmd2[DsizeL][2];
+  dataindexd2(level,out,din);
+  tmpmd2[0][0]=0;
+  for(k=1;k<=din[0][0];k++){
+    if(k!=ii){
+      js=din[k][0]; je=din[k][1];
+      appendd2(level,js,je,out,tmpmd2);
+    }else{
+      appendd2(level,1,nr,rep,tmpmd2);
+    }
+  }
+  nall=length2(tmpmd2);
+  out[0][0]=0;
+  appendd2(level,1,nall,tmpmd2,out);
+}
+
+void replacelistd3(int level,int ii,double out[][3],double rep[][3]){//180428
+  int din[DsizeS][2],k,js,je, nall=length3(out),nr=length3(rep);
+  double tmpmd3[DsizeL][3];
+  dataindexd3(level,out,din);
+  tmpmd3[0][0]=0;
+  for(k=1;k<=din[0][0];k++){
+    if(k!=ii){
+      js=din[k][0]; je=din[k][1];
+      appendd3(level,js,je,out,tmpmd3);
+    }else{
+      appendd3(level,1,nr,rep,tmpmd3);
+    }
+  }
+  nall=length3(tmpmd3);
+  out[0][0]=0;
+  appendd3(level,1,nall,tmpmd3,out);
+}
+
+void replacelistd6(int level,int ii,double out[][6],double rep[][6]){//180421
+  int din[DsizeS][2],k,js,je, nall=length6(out),nr=length6(rep);
+  double tmpmd6[DsizeL][6];
+  dataindexd6(level,out,din);
+  js=din[ii][0]; je=din[ii][1];
+  tmpmd6[0][0]=0;
+  for(k=1;k<=din[0][0];k++){
+    if(k!=ii){
+      js=din[k][0]; je=din[k][1];
+      appendd6(level,js,je,out,tmpmd6);
+    }else{
+      appendd6(level,1,nr,rep,tmpmd6);
+    }
+  }
+  nall=length6(tmpmd6);
+  out[0][0]=0;
+  appendd6(level,1,nall,tmpmd6,out);
+}
+
+void removelistd2(int level,int ii,double out[][2]){//180503
+  int din[DsizeS][2],k,js,je,nall=length2(out);
+  double tmpmd[DsizeL][2];
+  dataindexd2(level,out,din);
+  js=din[ii][0]; je=din[ii][1];
+  tmpmd[0][0]=0;
+  appendd2(0,je+1,nall,out,tmpmd);
+  out[0][0]=0;
+  if(js>1){
+    out[0][0]=js-2;
+  }
+  appendd2(0,1,length2(tmpmd),tmpmd,out);
+}
+
 void dispvec(int n,double dt[]){
  int i;
  for(i=0; i<n; i++){
@@ -433,6 +638,15 @@ void dispmatd6(int from, int upto, double mat[][6]){
   }
 }
 
+void dispmatd1all(double dt[]){
+ int i,n;
+ n=length1(dt);
+ for(i=1; i<=n; i++){
+   printf("%f ",dt[i]);
+ }
+ printf("\n");
+}
+
 void dispmatd2all(double mat[][2]){
   int i;
   if(mat[0][0]==0){printf("%s\n","no data");return;}
@@ -490,58 +704,68 @@ void dispmati1(int from, int upto, int mat[]){
   }
 }
 
+void writedata3(const char *fname, double out[][3]){
+  int i, nall;
+  FILE *fp;
+  nall=length3(out);
+  fp=fopen(fname,"w");
+  for(i=1; i<=nall; i++){
+    fprintf(fp,"%7.5f %7.5f %7.5f",out[i][0],out[i][1],out[i][2]);
+    if(i<nall){fprintf(fp,"\n");}
+  }
+  fclose(fp);
+}
+
 void readdata3(const char *fname, double out[][3]){
-  double data[DsizeL][4];
-  int i, j, ndata, nall;
+  int n,j;
   FILE *fp;
   fp=fopen(fname,"r");
   if(fp==NULL){
     printf("file not found\n");
     return;
   }
-  ndata=0;
-  i=1;
-  while( ! feof(fp) && i<DsizeL){
-    for(j=0;j<4;j++){
-      fscanf(fp,"%lf",&data[i][j]);
+  n=1;
+  while( ! feof(fp) && n<DsizeL){
+    for(j=0;j<3;j++){
+      fscanf(fp,"%lf",&out[n][j]);
     }
-    i++;
-    ndata++;
+    n++;
   }
-  ndata--;
+  n--;
   fclose(fp);
-  out[0][0]=0; nall=0;
-  for(i=1; i<=ndata; i++){
-    nall=add3(out, data[i][0], data[i][1], data[i][2]);
+  out[0][0]=n;
+}
+
+void writedata2(const char *fname, double out[][2]){
+  int i, nall;
+  FILE *fp;
+  nall=length2(out);
+  fp=fopen(fname,"w");
+  for(i=1; i<=nall; i++){
+    fprintf(fp,"%7.5f %7.5f",out[i][0],out[i][1]);
+    if(i<nall){fprintf(fp,"\n");}
   }
-  return;
+  fclose(fp);
 }
 
 void readdata2(const char *fname, double out[][2]){
-  double data[DsizeL][3];
-  int i, j, ndata, nall;
+  int n,j;
   FILE *fp;
   fp=fopen(fname,"r");
   if(fp==NULL){
     printf("file not found\n");
     return;
   }
-  ndata=0;
-  i=1;
-  while( ! feof(fp) && i<DsizeL){
-    for(j=0;j<3;j++){
-      fscanf(fp,"%lf",&data[i][j]);
+  n=1;
+  while( ! feof(fp) && n<DsizeL){
+    for(j=0;j<2;j++){
+      fscanf(fp,"%lf",&out[n][j]);
     }
-    i++;
-    ndata++;
+    n++;
   }
-  ndata--;
+  n--;
   fclose(fp);
-  out[0][0]=0; nall=0;
-  for(i=1; i<=ndata; i++){
-    nall=add2(out, data[i][0], data[i][1]);
-  }
-  return;
+  out[0][0]=n;
 }
 
 int output2(const char *var, const char *fname, int level, double out[][2]){
@@ -729,6 +953,7 @@ int connectseg(double pdata[][2], double eps, double out[][2]){
   int nall=length2(pdata), nout=0;
   double qd[DsizeM][2], ah[2], ao[2], p[2], q[2], tmp[2];
   out[0][0]=0;
+  if(length2(pdata)==0){return 0;} //180613
   nd=dataindexd2(0,pdata,din);
   qd[0][0]=0;
   nq=appendd2(0,din[1][0],din[1][1],pdata,qd);
@@ -780,6 +1005,63 @@ int connectseg(double pdata[][2], double eps, double out[][2]){
   return nout;
 }
 
+int connectseg3(double pdata[][3], double eps, double out[][3]){
+  int din[DsizeM][2],i,j,nd,nq,flg, st, en, ntmp;
+  int nall=length3(pdata), nout=0;
+  double qd[DsizeM][3], ah[3], ao[3], p[3], q[3], tmp[3];
+  out[0][0]=0;
+  if(length3(pdata)==0){return 0;} //180613
+  nd=dataindexd3(0,pdata,din);
+  qd[0][0]=0;
+  nq=appendd3(0,din[1][0],din[1][1],pdata,qd);
+  removei2(1,din);
+  while(nd>0){
+    pull3(1,qd,ah);
+    pull3(nq,qd,ao);
+	flg=0;
+    for(j=1; j<=nd; j++){
+	  st=din[j][0];
+      en=din[j][1];
+      pull3(st,pdata,p);
+      pull3(en,pdata,q);
+      if(dist(3,p,ao)<eps){
+        nq=appendd3(0,st+1,en,pdata,qd);
+        nd=removei2(j,din);
+        flg=1;
+        break;
+      }
+      if(dist(3,q,ao)<eps){
+        nq=appendd3(0,en-1,st,pdata,qd);
+        nd=removei2(j,din);
+        flg=1;
+        break;
+      }
+      if(dist(3,p,ah)<eps){
+        nq=prependd3(en,st+1,pdata,qd);
+        nd=removei2(j,din);
+        flg=1;
+        break;
+      }
+      if(dist(3,q,ah)<eps){
+        nq=prependd3(st,en-1,pdata,qd);
+        nd=removei2(j,din);
+        flg=1;
+        break;
+      }
+    }
+    if(flg==0){
+	  nout=appendd3(2,1,nq,qd,out);
+	  qd[0][0]=0;
+      nq=appendd3(0,din[1][0],din[1][1],pdata,qd);
+      nd=removei2(1,din);
+	}
+  }
+  if(nq>0){
+    nout=appendd3(2,1,nq,qd,out);
+  }
+  return nout;
+}
+
 void koutenseg(double a[2], double b[2], double c[2], double d[2],
           double eps, double eps2, double out[4]){
   double v[2], epss,eps0=Eps, sv2, p1, p2, q1,q2,m1, m2, M1, M2, t;
@@ -1227,26 +1509,26 @@ int dropnumlistcrv(double qd[][2], double eps1, int out[]){
     pd[0][0]=0; npd=0;
     npd=appendd2(0,se,en,qd,pd);
     ptL[0]=0; nptL=0;
-    nptL=appendi1(ptL,1);
+    nptL=appendi1(1,ptL);
     pull2(1,pd,p);
     for(k=2; k<=npd-1; k++){
       pull2(k,pd,tmp2d);
       if(dist(2,p,tmp2d)>eps1){
-        nptL=appendi1(ptL,k);
+        nptL=appendi1(k,ptL);
         pull2(k,pd,p);
       }
     }
     pull2(npd-1,pd,p);
     pull2(npd,pd,tmp2d);
     if(dist(2,p,tmp2d)>eps){  // eps -> eps1 ? 
-      nptL=appendi1(ptL,npd);
+      nptL=appendi1(npd,ptL);
     }
     if(nptL==1){ptL[0]=0; nptL=0;}
     if(nout>0){
-      nout=appendi1(out,Infint);
+      nout=appendi1(Infint,out);
     }
     for(i=1; i<=nptL; i++){
-      nout=appendi1(out,ptL[i]);
+      nout=appendi1(ptL[i],out);
     }
   }
   out[0]=nout;
@@ -1955,7 +2237,7 @@ void intersectpartseg(double crv1[][2],double crv2[][2],int ii, int jj,int num,d
 
 void collectsameseg(double ptdL[][6],double rL[][6]){
   double Eps00=pow(10,-8.0),tmp1d[2],tmp2d[2],tmp1d6[6],tmp2d6[6];
-  int ii,jj,kk,nr,ctr,numL[40],diffL[40];
+  int ii,jj,kk,istr,nr,numL[40],diffL[40];
   double s1,e1,s2,e2,dst,tmp1,tmp2,tmp1md6[40][6];
   rL[0][0]=0;
   if(length6(ptdL)==0){
@@ -1979,13 +2261,13 @@ void collectsameseg(double ptdL[][6],double rL[][6]){
   }else{
     s2=kk-Eps00; e2=s2+1+2*Eps00;
   }
-  ctr=1;numL[0]=0;
+  numL[0]=0;
   for(ii=1;ii<=length6(rL);ii++){
     pull6(ii,rL,tmp1d6);
     tmp1=tmp1d6[2]; tmp2=tmp1d6[3];
     if((tmp1>s1)&&(tmp1<e1)&&(tmp2>s2)&&(tmp2<e2)){
       addptd6(tmp1md6,tmp1d6);
-      numL[ctr]=ii;ctr++;
+      appendi1(ii,numL);
     }
   }
   dst=100;
@@ -2005,7 +2287,7 @@ void collectsameseg(double ptdL[][6],double rL[][6]){
   for(ii=0;ii<=length6(rL);ii++){
     diffL[ii]=0;
   }
-  for(ii=1;ii<=ctr;ii++){
+  for(ii=1;ii<=numL[0];ii++){
     kk=numL[ii];
     diffL[kk]=1;
   }
@@ -2018,12 +2300,12 @@ void collectsameseg(double ptdL[][6],double rL[][6]){
   }
 }
 
-void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][6]){
+void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int nbez,double out[][6]){
   double Eps00=pow(10,-8.0),Dist=10*Eps2, crv1[DsizeL][2],crv2[DsizeL][2];
-  double tmp1md6[DsizeS][6],tmp2md6[DsizeS][6];
+  double tmp1md6[DsizeS][6],tmp2md6[DsizeS][6],tmpmd1[DsizeS];
   double tmpd[2],tmp1d[2],tmp2d[2],tmpd6[6],tmp1d6[6];
-  double sx,sy,dst,ans[4];
-  int ii,jj,kk,nall,self,js,je;
+  double sx,sy,dst,ans1[4],ans2[4];
+  int ii,jj,kk,nall,self,js,je,din[DsizeS][2];
   crv1[0][0]=0;
   pull2(1,crv1s,tmpd);addptd2(crv1,tmpd);
   for(ii=2;ii<=length2(crv1s);ii++){
@@ -2063,7 +2345,7 @@ void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][
       js=ii+2; je=length2(crv2)-1;
     }
     for(jj=js;jj<=je;jj++){
-      intersectpartseg(crv1,crv2,ii,jj,num,tmpd6);
+      intersectpartseg(crv1,crv2,ii,jj,nbez,tmpd6);
       if(tmpd6[0]<Inf-Eps){
 		tmpd[0]=tmpd6[0];tmpd[1]=tmpd6[1];
         if(length6(tmp1md6)==0){
@@ -2085,41 +2367,48 @@ void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][
   tmp2md6[0][0]=0;
   while(length6(tmp1md6)>0){
 	collectsameseg(tmp1md6,tmp2md6);
-	dst=100;
-    for(ii=1;ii<=length6(tmp1md6);ii++){
-      pull6(ii,tmp1md6,tmpd6);
-      if(tmpd6[4]<dst){
-        dst=tmpd6[4];
-      }
+    appendd6(2,1,length6(tmp1md6),tmp1md6,out);
+	tmp1md6[0][0]=0;
+    for(ii=1;ii<=length6(tmp2md6);ii++){
+      pull6(ii,tmp2md6,tmpd6);addptd6(tmp1md6,tmpd6);
     }
-    sx=0; sy=0;nall=0;
-    for(ii=1;ii<=length6(tmp1md6);ii++){
-      pull6(ii,tmp1md6,tmpd6);
-      if(tmpd6[4]<dst+Eps00){
-        nall++;
-        sx=sx+tmpd6[0]; sy=sy+tmpd6[1];
+  }
+  dataindexd6(2,out,din);
+  for(ii=1;ii<=length2i(din);ii++){
+    tmp1md6[0][0]=0;
+    appendd6(2,din[ii][0],din[ii][1],out,tmp1md6);
+    if(length6(tmp1md6)==1){
+      replacelistd6(2,ii,out,tmp1md6);
+    }else{
+      dst=100;
+      for(jj=1;jj<=length6(tmp1md6);jj++){
+        pull6(jj,tmp1md6,tmpd6);
+        if(tmpd6[4]<dst){
+          dst=tmpd6[4];
+        }
+      }
+	  sx=0; sy=0;nall=0;
+      for(jj=1;jj<=length6(tmp1md6);jj++){
+        pull6(jj,tmp1md6,tmpd6);
+        if(tmpd6[4]<dst+Eps00){
+          nall++;
+          sx=sx+tmpd6[0]; sy=sy+tmpd6[1];
+        }
       }
       tmp1d[0]=sx/nall; tmp1d[1]=sy/nall;
-      tmp1d6[0]=tmp1d[0];tmp1d6[1]=tmp1d[1];
-      nearestptpt(tmp1d,crv1,ans);
-      tmp1d6[2]=ans[2];
-      nearestptpt(tmp1d,crv2,ans);
-      tmp1d6[3]=ans[2];
-      tmp1d6[4]=dst;
-      tmp1d6[5]=tmpd6[5];
-      addptd6(out,tmp1d6);
-    }
-    tmp1md6[0][0]=0;
-    for(ii=1;ii<=length6(tmp2md6);ii++){
-      pull6(ii,tmp2md6,tmpd6);addptd6(tmp1md6,tmpd6);
+      nearestptpt(tmp1d,crv1,ans1);
+      nearestptpt(tmp1d,crv2,ans2);
+      tmp2md6[0][0]=0;
+      push6(sx/nall,sy/nall,ans1[2],ans2[2],dst,tmp1md6[1][5],tmp2md6);
+      replacelistd6(2,ii,out,tmp2md6);
     }
   }
 }
 
-void intersectcurves(double crv1[][2],double crv2[][2],int num,double out[][2]){
+void intersectcurves(double crv1[][2],double crv2[][2],int nbez,double out[][2]){
   double out6[DsizeS][6],tmpd[2],tmpd6[6];
   int jj;
-  intersectcurvesPp(crv1,crv2,num,out6);
+  intersectcurvesPp(crv1,crv2,nbez,out6);
   out[0][0]=0;
   for(jj=1;jj<=length6(out6);jj++){
     pull6(jj,out6,tmpd6);push2(tmpd6[0],tmpd6[1],out);
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h
index 7e83a2b9e83..7b22ddfba83 100644
--- a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h
@@ -12,6 +12,7 @@ int addptd3(double gdata[][3], double pt[3]);
 int addptd4(double gdata[][4], double pt[4]);
 int addptd5(double gdata[][5], double pt[5]);
 int addptd6(double gdata[][6], double pt[6]);
+int add1(double gdata[], double x);
 int add2(double gdata[][2], double x, double y);
 int add3(double gdata[][3], double x, double y, double z);
 int push2(double x, double y, double ptL[][2]);
@@ -26,18 +27,29 @@ void pull5(int num, double ptL[][5],double pt[5]);
 void pull6(int num, double ptL[][6],double pt[6]);
 int appendd2(int level, int from,int upto, double gdata[][2],double mat[][2]);
 int appendd3(int level, int from,int upto, double gdata[][3],double mat[][3]);
+int appendd6(int level, int from,int upto, double gdata[][6],double mat[][6]);
 int prependd2(int from,int upto, double gdata[][2],double mat[][2]);
-int appendd1(double numvec[], double num);
-int appendi1(int numvec[],int num);
+int prependd3(int from,int upto, double gdata[][3],double mat[][3]);
+int appendd1(double num,double numvec[]);
+int prependd1(double num, double numvec[]);
+int appendi1(int num,int numvec[]);
 int prependi1(int num, int numvec[]);
 int removei2(int nrmv,int mat[][2]);
+void removed3(int nrmv,double mat[][3]);
 int dataindexi1(int out[],int din[][2]);
 int dataindexd1(double out[],int din[][2]);
 int dataindexd2(int level,double out[][2],int din[][2]);
 int dataindexd3(int level, double out[][3],int din[][2]);
 int dataindexd4(int level, double out[][4],int din[][2]);
+int dataindexd6(int level, double out[][6],int din[][2]);
+void replacelistd1(int ii,double out[],double rep[]);
+void replacelistd2(int level,int jj,double out[][2],double rep[][2]);
+void replacelistd3(int level,int jj,double out[][3],double rep[][3]);
+void replacelistd6(int level,int ii,double out[][6],double rep[][6]);
+void removelistd2(int level,int ii,double out[][2]);
 void dispvec(int n,double dt[]);
 void dispveci(int n,int dt[]);
+void dispmatd1all(double dt[]);
 void dispmatd2(int from, int upto, double mat[][2]);
 void dispmatd3(int from, int upto, double mat[][3]);
 void dispmatd4(int from, int upto, double mat[][4]);
@@ -50,7 +62,9 @@ void dispmatd5all(double mat[][5]);
 void dispmatd6all(double mat[][6]);
 void dispmati2(int from, int upto, int mat[][2]);
 void dispmati1(int from, int upto, int mat[]);
+void writedata3(const char *fname, double out[][3]);
 void readdata3(const char *fname, double out[][3]);
+void writedata2(const char *fname, double out[][2]);
 void readdata2(const char *fname, double out[][2]);
 int output2(const char *var, const char *fname, int level, double out[][2]);
 void simplesort(double number[]);
@@ -64,6 +78,7 @@ void pointoncurve(double t, double gdata[][2], double eps, double pt[]);
 int partcrv(double a,double b,double pkL[][2], double pL[][2]);
 int partcrv3(double a,double b,double pkL[][3], double pL[][3]);
 int connectseg(double pdata[][2], double eps, double out[][2]);
+int connectseg3(double pdata[][3], double eps, double out[][3]);
 void koutenseg(double a[2], double b[2], double c[2], double d[2],
           double eps, double eps2, double out[4]);
 double paramoncurve(double p[2], int n, double ptL[][2]);
@@ -90,8 +105,8 @@ void intersectseg(double seg1[][2],double seg2[][2], double out[5]);
 int osplineseg(double ptlist[5][2],int num,double out[][2]);
 void intersectpartseg(double crv1[][2],double crv2[][2],int ii, int jj,int num,double out[6]);
 void collectsameseg(double ptdL[][6],double rL[][6]);
-void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][6]);
-void intersectcurves(double crv1[][2],double crv2[][2],int num,double out[][2]);
+void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int nbez,double out[][6]);
+void intersectcurves(double crv1[][2],double crv2[][2],int nbez,double out[][2]);
 
 const double Inf=99999;
 const int Infint=30000;
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketlibClog.txt b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketlibClog.txt
new file mode 100644
index 00000000000..414cc1e0346
--- /dev/null
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketlibClog.txt
@@ -0,0 +1,18 @@
+180613
+  connectseg,connectseg3 changed ( for null data )
+  output3h debugged  (case where din1[0][0]==1)
+180609
+  crv3onsfparadata debugged ( case where curve has multi parts)
+  connectseg3 added
+180601
+  wireparadata changed (for cmdflg=1)
+  sfcutparadata changed ( for cmdflg=1)
+180531
+  outputend added
+  output3,output3h changed
+  readoutdata3 changed  ( varname added )
+  intersectcrvsf changed (ws added)
+180530
+  readoutdata3 added ( hidden data not included )
+180523
+  crvsfparadata debugged (npar removed)
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h
index d9c0fbdc18c..7f5d9923c9a 100644
--- a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h
@@ -1,19 +1,24 @@
-// 17.06.16 crvsfparadata changed ( num of devision )
+// 180430 
+//    crv2sfparadata,wireparadata added
+//    output3h,output3 chanded       
+// 180426 readdataC,writedataC changed
+// 180420 newly rewrited
+// 170616 crvsfparadata changed ( num of devision )
 
-int output3(const char *var, const char *fname, int level, double out[][3]){
-  int nall=floor(out[0][0]+0.5), ndin, din[DsizeM][2],i,j,ctr;
+int output3(short head, const char *wa,const char *var, const char *fname,double out[][3]){
+  int din[DsizeS][2],i,j,ctr;
   double tmpd[3];
   FILE *fp;
-  fp=fopen(fname,"w");
+  fp=fopen(fname,wa);
   if (fp == NULL) { 
     printf("cannot open\n");  
     return -1;
-//    exit(1); 
   }
-  ndin=dataindexd3(level, out,din);
-  fprintf(fp,"%s//\n",var); //17.05.29
-//  fprintf(fp,"%s3d//\n",var2);
-  for(i=1; i<=ndin; i++){
+  dataindexd3(2, out,din);
+  if(head==1){
+    fprintf(fp,"%s//\n",var);
+  }
+  for(i=1; i<=length2i(din); i++){
     fprintf(fp,"start//\n");
     fprintf(fp,"[");
     ctr=0;
@@ -22,32 +27,31 @@ int output3(const char *var, const char *fname, int level, double out[][3]){
       fprintf(fp,"[ %7.5f, %7.5f, %7.5f]",tmpd[0],tmpd[1],tmpd[2]);
       ctr++;
       if(ctr==5){
-        fprintf(fp,"]//\n[");
+        fprintf(fp,"]//\n");
         ctr=0;
+        if(j<din[i][1]){
+          fprintf(fp,"[");
+        }
       }else{
         if(j<din[i][1]){
           fprintf(fp,",");
+        }else{
+          fprintf(fp,"]//\n");
         }
       }
     }
-    fprintf(fp,"]//\n");
-    if(i<ndin){
-      fprintf(fp,"end//\n");
-    }
-    else{
-      fprintf(fp,"end////\n");
-    }
+    fprintf(fp,"end//\n");
   }
   fclose(fp);
   return 0;
 }
 
-int output3h(const char *var, const char *varh, const char *fname, double out[][3]){
+int output3h(const char *wa,const char *var, const char *varh, const char *fname, double out[][3]){
   // 17.06.09
   double out1[DsizeLL][3], tmpd[3];
   int din1[DsizeM][2],din2[DsizeM][2], i, j,ctr;
   FILE *fp;
-  fp=fopen(fname,"w");
+  fp=fopen(fname,wa);
   if (fp == NULL) { 
     printf("cannot open\n");  
     return -1;
@@ -81,7 +85,9 @@ int output3h(const char *var, const char *varh, const char *fname, double out[][
   }
   out1[0][0]=0;
   fprintf(fp,"%s//\n",varh); //17.06.16from
-  appendd3(0, din1[2][0],din1[2][1],out,out1);
+  if(din1[0][0]>1){ //180613from
+    appendd3(0, din1[2][0],din1[2][1],out,out1);
+  }//180613to
   dataindexd3(2, out1,din2);
   for(i=1; i<=length2i(din2); i++){
     fprintf(fp,"start//\n");
@@ -102,52 +108,54 @@ int output3h(const char *var, const char *varh, const char *fname, double out[][
       }
     }
     fprintf(fp,"]//\n");
-//    if(i<ndin2){
-      fprintf(fp,"end//\n");
-//    }
-//    else{
-//      fprintf(fp,"end////\n");
-//    }
+    fprintf(fp,"end//\n");
   }
-  fprintf(fp,"//\n");
+//  fprintf(fp,"//\n"); //180531
   fclose(fp);
   return 0;
 }
 
+void outputend(const char *fname){
+  FILE *fp;
+  fp=fopen(fname,"a");
+  fprintf(fp,"//\n");
+  fclose(fp);
+}
+
 int writedataC(const char *fname, double out[][3]){
   int i, nall;
   FILE *fp;
   nall=length3(out);
+  printf("123 %d\n",nall);
   fp=fopen(fname,"w");
   for(i=1; i<=nall; i++){
-    fprintf(fp,"%7.5f %7.5f %7.5f\n",out[i][0],out[i][1],out[i][2]);
+    fprintf(fp,"%7.5f %7.5f %7.5f %6s",out[i][0],out[i][1],out[i][2],"-99999");
+    if(i<nall){
+      fprintf(fp,"\n");
+    };
   }
-  fprintf(fp,"%7.5f %7.5f %7.5f\n",Inf,Inf,Inf);
   fclose(fp);
   return nall;
 }
 
 int readdataC(const char *fname, double data[][3]){
   int i, j, ndata=0;
+  double tmpd4[4];
   FILE *fp;
   fp=fopen(fname,"r");
   if(fp==NULL){
     printf("file not found\n");
     return 1;
   }
+  data[0][0]=0;
   while( ! feof(fp) && ndata<DsizeL){
-    ndata++;
-    for(j=0;j<=2;j++){
-      fscanf(fp,"%lf",&data[ndata][j]);
-    }
-    if(data[ndata][1]==Inf){
-      ndata--;
-      break;
+    for(j=0;j<=3;j++){
+      fscanf(fp,"%lf",&tmpd4[j]);
     }
+    push3(tmpd4[0],tmpd4[1],tmpd4[2],data);
   }
   fclose(fp);
-   data[0][0]=ndata;
-  return ndata;
+  return length3(data);
 }
 
 void parapt(double pt3[3],double pt2[2]){
@@ -160,6 +168,12 @@ void parapt(double pt3[3],double pt2[2]){
   }
 }
 
+double zparapt(double p[3]){
+  double x=p[0], y=p[1], z=p[2], out;
+  out=x*cos(PHI)*sin(THETA)+y*sin(PHI)*sin(THETA)+z*cos(THETA);
+  return out;
+}
+
 int projpara(double p3[][3], double p2[][2]){
   double pt3[3], pt2[2];
   int i, nall, nall2;
@@ -198,67 +212,17 @@ double invparapt(double ph, double fh[][2], double fk[][3],double out[3]){
   return n+s0;
 }
 
-int surfcurve(short ch,int level, double crv[][2], double pdata[][3]){
-  double p2[2], pt[3];
-  int i, n, nall;
-  n=length2(crv);
+void surfcurve(short ch,double crv[][2], double pdata[][3]){
+  double pd2[2], pt[3];
+  int i;
   pdata[0][0]=0;
-  for(i=1; i<=n; i++){
-    pull2(i,crv,p2);
-    surffun(ch,p2[0],p2[1],pt);
-    nall=addptd3(pdata,pt);
+  for(i=1; i<=length2(crv); i++){
+    pull2(i,crv,pd2);
+    surffun(ch,pd2[0],pd2[1],pt);
+    addptd3(pdata,pt);
   }
-  return nall;
 }
 
-int xyzax3data(double x0, double x1, double y0, double y1, double z0, double z1, 
-                           double out[][3]){
-  int nall;
-  double pd3[10][3];
-  out[0][0]=0; nall=0;
-  pd3[0][0]=2;
-  pd3[1][0]=x0; pd3[1][1]=0; pd3[1][2]=0;
-  pd3[2][0]=x1; pd3[2][1]=0; pd3[2][2]=0;
-  nall=appendd3(2,1,2,pd3,out);
-  pd3[1][0]=0; pd3[1][1]=y0; pd3[1][2]=0;
-  pd3[2][0]=0; pd3[2][1]=y1; pd3[2][2]=0;
-  nall=appendd3(2,1,2,pd3,out);
-  pd3[1][0]=0; pd3[1][1]=0; pd3[1][2]=z0;
-  pd3[2][0]=0; pd3[2][1]=0; pd3[2][2]=z1;
-  nall=appendd3(2,1,2,pd3,out);
-  return nall;
-}
-
-void evlptablepara(short ch){
-  int i, j;
-  double u,v,u1,u2,v1,v2,pt1[3],pt2[3],pt2d1[2],pt2d2[2],dxu,dyu,dxv,dyv;
-  double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
-  double eps=0.00001;
-  for (j= 0; j < Ndv+1; j++) {
-    v=Vrng[0]+j*dv;
-    v1=v-eps/2; v2=v+eps/2;
-    for (i = 0; i < Mdv+1; i++) {
-      u=Urng[0]+i*du;
-      u1=u-eps/2; u2=u+eps/2;
-      surffun(ch,u1,v,pt1); parapt(pt1,pt2d1);
-	  surffun(ch,u2,v,pt2); parapt(pt2,pt2d2);	  
-      dxu=(pt2d2[0]-pt2d1[0])/eps;
-      dyu=(pt2d2[1]-pt2d1[1])/eps;
-      surffun(ch,u,v1,pt1); parapt(pt1,pt2d1);
-      surffun(ch,u,v2,pt2); parapt(pt2,pt2d2);
-      dxv=(pt2d2[0]-pt2d1[0])/eps;
-      dyv=(pt2d2[1]-pt2d1[1])/eps;
-      Zval[j][i]=dxu*dyv-dxv*dyu;     
-    }
-  }
-  for (i = 0; i < Mdv+1; i++) {
-    Xval[i]=Urng[0]+i*du;
-  }
-  for (j= 0; j < Ndv+1; j++) {
-    Yval[j]=Vrng[0]+j*dv;
-  }
-};
-
 double evlpfun(short ch, double u, double v){
   double u1,u2,v1,v2,pt1[3],pt2[3],pt2d1[2],pt2d2[2],dxu,dyu,dxv,dyv;
   double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
@@ -285,7 +249,7 @@ int envelopedata(short ch, double out[][2]){
   double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
   double pl[5][2], vl[5], ql[11][2], out2[DsizeM][2];
   double du=(Urng[1]-Urng[0])/Mdv;
-  double dv=(Vrng[1]-Vrng[0])/Mdv;
+  double dv=(Vrng[1]-Vrng[0])/Ndv; //180427
   ctr=0;
   for (j = 0; j < Ndv; j++) {
     vval1=Vrng[0]+j*dv;
@@ -351,67 +315,6 @@ int envelopedata(short ch, double out[][2]){
   return nall;
 }
 
-/*
-
-int implicitplot(double out[][2]){
-  int i, j, k, ctrq=0,ctr, nall;
-  double pl[5][2], vl[5], ql[11][2], out2[DsizeM][2];
-  ctr=0;
-  for (j = 0; j < Ndv; j++) {
-    for (i = 0; i < Mdv; i++) {
-	  pl[0][0]=Xval[i]; pl[0][1]=Yval[j]; vl[0]=Zval[j][i];
-      pl[1][0]=Xval[i+1]; pl[1][1]=Yval[j]; vl[1]=Zval[j][i+1];
-      pl[2][0]=Xval[i+1]; pl[2][1]=Yval[j+1]; vl[2]=Zval[j+1][i+1];
-      pl[3][0]=Xval[i]; pl[3][1]=Yval[j+1]; vl[3]=Zval[j+1][i];
-	  pl[4][0]=Xval[i]; pl[4][1]=Yval[j]; vl[4]=Zval[j][i];
-      ctrq=0;
-      for (k = 0; k < 4; k++) {
-        if(fabs(vl[k])<=Eps){
-          ql[ctrq][0]=pl[k][0];
-          ql[ctrq][1]=pl[k][1];
-          ctrq++;
-        }
-        else if(vl[k]>Eps){
-          if(vl[k+1]< -Eps){
-            ql[ctrq][0]=1/(vl[k]-vl[k+1])*(-vl[k+1]*pl[k][0]+vl[k]*pl[k+1][0]);
-            ql[ctrq][1]=1/(vl[k]-vl[k+1])*(-vl[k+1]*pl[k][1]+vl[k]*pl[k+1][1]);
-            ctrq++;
-          }
-        }
-        else{
-          if(vl[k+1]>Eps){
-            ql[ctrq][0]=1/(-vl[k]+vl[k+1])*(vl[k+1]*pl[k][0]-vl[k]*pl[k+1][0]);
-            ql[ctrq][1]=1/(-vl[k]+vl[k+1])*(vl[k+1]*pl[k][1]-vl[k]*pl[k+1][1]);
-            ctrq++;
-          }
-        }
-      }
-      if(ctrq==2){
-        if(ctr>0){
-          ctr++;
-          out[ctr][0]=Inf;
-          out[ctr][1]=2;  // 2 <- 0 ?
-        }
-        ctr++;
-        out[ctr][0]=ql[0][0];
-        out[ctr][1]=ql[0][1];
-        ctr++;
-        out[ctr][0]=ql[1][0];
-        out[ctr][1]=ql[1][1];
-      }
-    }
-  }
-  out[0][0]=ctr;
-  if(ctr>0){
-    nall=connectseg(out, Eps, out2);
-    out[0][0]=0;
-    nall=appendd2(2,1,nall,out2,out);
-  }
-  return nall;
-}
-
-*/
-
 int cuspsplitpara(short ch,double pdata[][2], double outkL[][3]){
   int i, j, ng, n, nd, ndd, ndrop, is, ncusp=0, ncusppt,st, en, npk,nph, cflg;
   int cusp[DsizeM], cuspflg;
@@ -508,24 +411,24 @@ int cuspsplitpara(short ch,double pdata[][2], double outkL[][3]){
           ntmp=trunc((i+j)*0.5);
           i=j;
           if(cusp[0]==0){
-            ncusp=appendi1(cusp,ntmp);
+            ncusp=appendi1(ntmp,cusp);
           }
           else{
-            ncusp=appendi1(cusp,ntmp);
+            ncusp=appendi1(ntmp,cusp);
           }
         }
       }
     }
     if(cflg==0){
       ncusp=prependi1(1,cusp);
-      ncusp=appendi1(cusp,nph);
+      ncusp=appendi1(nph,cusp);
     }
     else if(ncusp==0){
        ncusp=prependi1(1,cusp);
-       ncusp=appendi1(cusp,nph);
+       ncusp=appendi1(nph,cusp);
     }
     else if(cusp[1]==1){
-      ncusp=appendi1(cusp,nph);
+      ncusp=appendi1(nph,cusp);
     }
     else{
       ntmp=cusp[1];
@@ -543,7 +446,7 @@ int cuspsplitpara(short ch,double pdata[][2], double outkL[][3]){
       for(i=1; i<=cusp[0]; i++){
         cusp[i]=cusp[i]-ntmp+1;
       }
-      ncusp=appendi1(cusp,nph);
+      ncusp=appendi1(nph,cusp);
     }
     if(ncusp==2){
       pull2(nph,pth,tmp);
@@ -579,7 +482,7 @@ int cuspsplitpara(short ch,double pdata[][2], double outkL[][3]){
     ntmp=appendd3(0,din[i][0],din[i][1],outkL,tmp3d);
 	tmpnum[0]=0; ntmpnum=0;
 	for(j=dind[i][0]; j<=dind[i][1]; j++){
-      ntmpnum=appendi1(tmpnum,drop[j]);
+      ntmpnum=appendi1(drop[j],tmpnum);
     }
     tmp3d1[0][0]=0; ntmp1=0;
     for(j=1; j<=ntmpnum; j++){
@@ -639,32 +542,32 @@ int makexybdy(short ch,double ehL[][2]){
 }
 
 void partitionseg(double fig[][2],double gL[][2], int ns, double parL[]){
-  int ii,jj,kk,n,din[DsizeM][2],din2[DsizeM][2];
+  int ii,jj,kk,n,din[DsizeM][2];
   double eps0=pow(10,-5.0),tmpd4[4], tmpd6[6],tmp,tmp2;
   double p[2], q[2], pd3[3], ans[4], tmp1d1[DsizeM], tmp2d1[DsizeM];
   double other[DsizeM], g[DsizeM][2], kouL[DsizeS][4],kouL6[DsizeS][6];
-  if(ns>1){
+  if(ns>0){
     dataindexd1(Otherpartition,din);
     other[0]=0;
     for(ii=din[ns][0]; ii<=din[ns][1]; ii++){
-      appendd1(other,Otherpartition[ii]);
+      appendd1(Otherpartition[ii],other);
     }
   }else{
     other[0]=0;
   }
   parL[0]=0;
-  appendd1(parL,1);
-  appendd1(parL,length2(fig));
+  appendd1(1,parL);
+  appendd1(length2(fig),parL);
   if(length1(other)>0){
     for(ii=1; ii<=length1(other); ii++){
-      appendd1(parL,other[ii]);
+      appendd1(other[ii],parL);
     }
   }
   dataindexd2(2,gL,din);
-  for(n=ns; n<=din[0][0]; n++){
+  for(n=fmax(ns,1); n<=din[0][0]; n++){
     g[0][0]=0; 
     appendd2(2,din[n][0],din[n][1],gL,g);
-    intersectcurvesPp(fig,g,20,kouL6);
+	intersectcurvesPp(fig,g,20,kouL6);
     kouL[0][0]=0;
     for(ii=1;ii<=length6(kouL6);ii++){
       pull6(ii,kouL6,tmpd6);
@@ -675,58 +578,42 @@ void partitionseg(double fig[][2],double gL[][2], int ns, double parL[]){
       pull4(jj,kouL,tmpd4);
       tmp=tmpd4[2];
       if((tmp>1+Eps)&&(tmp<length2(fig)-Eps)){
-        appendd1(tmp1d1,tmp);
+        appendd1(tmp,tmp1d1);
       }
       tmp=tmpd4[3];
       if((tmp>1+Eps)&&(tmp<length2(g)-Eps)){
-        appendd1(tmp2d1,tmp);
+        appendd1(tmp,tmp2d1);
       }
     }
     for(jj=1;jj<=length1(tmp1d1);jj++){
-      appendd1(parL,tmp1d1[jj]);
+      appendd1(tmp1d1[jj],parL);
     }
     if((ns>0)&&(n>ns)&&(length1(tmp2d1)>0)){
       dataindexd1(Otherpartition,din);
       tmp1d1[0]=0;
-      for(ii=1;ii<=length2i(din);ii++){
-        if(ii==n){
-          for(kk=1;kk<=length1(tmp2d1);kk++){
-            appendd1(tmp1d1,tmp2d1[kk]);
-          }
-        }
-        for(jj=din[ii][0]; jj<=din[ii][1];jj++){
-           appendd1(tmp1d1,Otherpartition[jj]);
-        }
-        if(ii<length2i(din)){
-          appendd1(tmp1d1,Inf);
-        }
-      }
-    }
-    dataindexd1(tmp1d1,din);
-    Otherpartition[0]=0;
-    for(ii=1;ii<=length2i(din);ii++){
-      for(jj=din[ii][0];jj<=din[ii][1];jj++){
-        appendd1(Otherpartition,tmp1d1[jj]);
+      for(jj=din[n][0];jj<=din[n][1];jj++){
+        appendd1(Otherpartition[jj],tmp1d1);
       }
+      replacelistd1(n,Otherpartition,tmp2d1);
     }
   }
   simplesort(parL);
   tmp1d1[0]=0;
   for(ii=1;ii<=length1(parL);ii++){
-    appendd1(tmp1d1,parL[ii]);
+    appendd1(parL[ii],tmp1d1);
   }
   parL[0]=0;
   tmp2=-1;
   for(jj=1;jj<=length1(tmp1d1);jj++){
     tmp=tmp1d1[jj];
     if(fabs(tmp-tmp2)>Eps1){
-      appendd1(parL,tmp);
+      appendd1(tmp,parL);
       tmp2=tmp;
     }
   }
 }
 
-double funpthiddenQ(short ch,double u, double v, double pa[3]){
+double funpthiddenQ(short ch,double u, double v,double pa[3]){
   double Vec[3]={sin(THETA)*cos(PHI),sin(THETA)*sin(PHI),cos(THETA)};
   double eps0=pow(10,-4.0)*(XMAX-XMIN)/10, pt[3], out;
   surffun(ch,u, v, pt);
@@ -739,248 +626,120 @@ double funpthiddenQ(short ch,double u, double v, double pa[3]){
   return out;
 }
 
-void makevaltable(double pta[3]){
-  double eps0=pow(10,-6.0), u, v, du, dv, p[3];
-  int i, j, nx=0, ny=0; ;
-  du=(Urng[1]-Urng[0])/Mdv;
-  dv=(Vrng[1]-Vrng[0])/Ndv;
-  for(j=0; j<=Ndv; j++){
-    v=Vrng[0]+j*dv;
-    for(i=0; i<=Mdv; i++){
-      u=Urng[0]+i*du;
-      Zval[j][i]=funpthiddenQ(Ch,u, v, pta);
-    }
-  }
-  for(i=0; i<=Mdv; i++){
-    Xval[i]=Urng[0]+i*du;
-  }
-  for(j=0; j<=Ndv; j++){
-    Yval[j]=Vrng[0]+j*dv;
-  }
-}
-
-int pthiddenQ(short ch,double pta[3], int uveq, double eps[2]){
-  int i, j, k, out, nqL, nsL, nptL, flg, nxt;
+double pthiddenQ(short ch,double pta[3], int uveq){
+  int i, j, k;
   double Vec[3]={sin(THETA)*cos(PHI),sin(THETA)*sin(PHI),cos(THETA)};
-  double eps0, epspt, vec1[3], po[3]={0,0,0}, tmpd1, tmp1d1, tmp2d1, m1, m2;
-  double tmp3d1, tmp4d1, tmp5d1, tmp6d1, tmpd3[3], v1, v2, v3;
-  double tmpd2[2], p1[2], p2[2], ptuv[2], puv1[2], puv2[2];
-  double  p1d3[3], p2d3[3], ptd3[3];
-  double ptL[DsizeM][6], sL[DsizeM][7], qL[DsizeM][2], pL[5][2], vL[5];
+  double v1,v2,v3,out=-Inf;
+  double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
+  double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
+  double pl[5][2], vl[5], ql[11][2],p1[2],p2[2],m1,m2;
+  double ptuv[2], puv1[2], puv2[2],p1d3[3], p2d3[3], ptd3[3];
+  double tmpd1, tmpd2[2],tmpd3[3],tmp1d3[3];
 //     Out=1 : hidden
 //     SL : List of segments near PtA+Vec
-  eps0=pow(10,-4.0)*(XMAX-XMIN)/10.0;
-  epspt=eps[0]*(XMAX-XMIN)/10.0;
-  out=0;
-  po[0]=0; po[1]=0; po[2]=0;
-  tmpd1=1/dist(3,Vec,po);
-  vec1[0]=tmpd1*Vec[0]; 
-  vec1[1]=tmpd1*Vec[1]; 
-  vec1[2]=tmpd1*Vec[2]; 
-  nptL=0;
-  makevaltable(pta);
-  nsL=0;
   for(j=0; j<=Ndv-1; j++){
+    vval1=Vrng[0]+j*dv;
+    vval2=Vrng[0]+(j+1)*dv;
+    uval1=Urng[0];
+    eval11=funpthiddenQ(ch,uval1,vval1,pta);
+    eval12=funpthiddenQ(ch,uval1,vval2,pta);
     for(i=0; i<=Mdv-1; i++){
-      pL[0][0]=Xval[i]; pL[0][1]=Yval[j]; vL[0]=Zval[j][i];
-      pL[1][0]=Xval[i+1]; pL[1][1]=Yval[j]; vL[1]=Zval[j][i+1];
-      pL[2][0]=Xval[i+1]; pL[2][1]=Yval[j+1]; vL[2]=Zval[j+1][i+1];
-      pL[3][0]=Xval[i]; pL[3][1]=Yval[j+1]; vL[3]=Zval[j+1][i];
-      pL[4][0]=Xval[i]; pL[4][1]=Yval[j]; vL[4]=Zval[j][i];
-      qL[0][0]=0; nqL=0;
+      uval2=Urng[0]+(i+1)*du;
+      eval21=funpthiddenQ(ch,uval2,vval1,pta);
+      eval22=funpthiddenQ(ch,uval2,vval2,pta);
+	  pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
+      pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
+      pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
+      pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
+	  pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
+      ql[0][0]=0;
       for(k=0; k<=3; k++){
-        p1[0]=pL[k][0]; p1[1]=pL[k][1]; 
-        p2[0]=pL[k+1][0]; p2[1]=pL[k+1][1]; 
-        m1=vL[k]; m2=vL[k+1];
-        if((m1>=eps0)&&(m2>=eps0)){
+        pull2(k,pl,p1);
+        pull2(k+1,pl,p2);
+        m1=vl[k]; m2=vl[k+1];
+        if(fabs(m1)<Eps){
+          addptd2(ql, p1);
           continue;
         }
-        if((m1<=-eps0)&&(m2<=-eps0)){
+        if(fabs(m2)<Eps){
           continue;
         }
-        if(fabs(m1)<eps0){
-          tmpd2[0]=p1[0]; tmpd2[1]=p1[1];
-        }
-        else if(fabs(m2)<eps0){
+        if((m1>0)&&(m2>0)){
           continue;
         }
-        else{
-          tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
-          tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
+        if((m1< 0)&&(m2< 0)){
+          continue;
         }
-        nqL++;
-        qL[nqL][0]=tmpd2[0]; qL[nqL][1]=tmpd2[1];
+        tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
+        tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
+        addptd2(ql, tmpd2);
       }
-      if(nqL==2){
-        puv1[0]=qL[1][0]; puv1[1]=qL[1][1]; 
-        puv2[0]=qL[2][0]; puv2[1]=qL[2][1];  
+      uval1=uval2;
+      eval11=eval21;
+      eval12=eval22;
+      if(length2(ql)==2){
+        puv1[0]=ql[1][0]; puv1[1]=ql[1][1]; 
+        puv2[0]=ql[2][0]; puv2[1]=ql[2][1];  
         surffun(ch,puv1[0],puv1[1],p1d3);
         surffun(ch,puv2[0],puv2[1],p2d3);
-        if(dist(2,p1d3,p2d3)<eps0){
-          continue;
-        }
-        tmp1d1=dist(3,p1d3,pta);
-        tmpd3[0]=p1d3[0]-pta[0];
-        tmpd3[1]=p1d3[1]-pta[1];
-        tmpd3[2]=p1d3[2]-pta[2]; 
-        tmp2d1=dotprod(3,tmpd3,vec1);
-        tmp3d1=sqrt(fabs(pow(tmp1d1,2.0)-pow(tmp2d1,2.0)));
-        tmp4d1=dist(3,p2d3, pta);
-        tmpd3[0]=p2d3[0]-pta[0];
-        tmpd3[1]=p2d3[1]-pta[1];
-        tmpd3[2]=p2d3[2]-pta[2]; 
-        tmp5d1=dotprod(3,tmpd3,vec1);
-        tmp6d1=sqrt(fabs(pow(tmp4d1,2.0)-pow(tmp5d1,2.0)));
-        if(tmp2d1<tmp5d1){
-          tmpd3[0]=p1d3[0]; tmpd3[1]=p1d3[1]; tmpd3[2]=p1d3[2];
-          p1d3[0]=p2d3[0]; p1d3[1]=p2d3[1]; p1d3[2]=p2d3[2];
-          p2d3[0]=tmpd3[0]; p2d3[1]=tmpd3[1]; p2d3[2]=tmpd3[2];
-          tmpd1=tmp2d1; tmp2d1=tmp5d1; tmp5d1=tmpd1;
-          tmpd1=tmp3d1; tmp3d1=tmp6d1; tmp6d1=tmpd1;
-        }
-        if((tmp3d1<eps[1]*epspt)&&(tmp6d1<eps[1]*epspt)){
-          if(tmp2d1>-eps0){
-            if(tmp5d1<-eps0){
-              nsL++;
-              sL[nsL][0]=p1d3[0]; sL[nsL][1]=p1d3[1]; sL[nsL][2]=p1d3[2];
-              sL[nsL][3]=p2d3[0]; sL[nsL][4]=p2d3[1]; sL[nsL][5]=p2d3[2];
-              sL[nsL][6]=-1;
-            }
-            else{
-              nsL++;
-              sL[nsL][0]=p1d3[0]; sL[nsL][1]=p1d3[1]; sL[nsL][2]=p1d3[2];
-              sL[nsL][3]=p2d3[0]; sL[nsL][4]=p2d3[1]; sL[nsL][5]=p2d3[2];
-              sL[nsL][6]=1;
-            }
-          }
-        }
-        v1=vec1[0]; v2=vec1[1]; v3=vec1[2];
-        if(fabs(v1)>eps0){
+        v1=Vec[0]; v2=Vec[1]; v3=Vec[2];
+        if(fabs(v1)>Eps){
           m1=pta[2]+v3/v1*(p1d3[0]-pta[0])-p1d3[2];
           m2=pta[2]+v3/v1*(p2d3[0]-pta[0])-p2d3[2];
-        }
-        else if(fabs(v2)>eps0){
+        }else if(fabs(v2)>Eps){
           m1=pta[2]+v3/v2*(p1d3[1]-pta[1])-p1d3[2];
           m2=pta[2]+v3/v2*(p2d3[1]-pta[1])-p2d3[2];
-        }
-        else{
+        }else{
           m1=pta[1]-p1d3[1];
-          m2=pta[1]-p2d3[1];
+          m2=pta[1]-p2d3[2];
         }
-        if(m1*m2>=0){
-          if(((m1>eps0)&&(m2>eps0))||((m1<-eps0)&&(m2<-eps0))){
+        if(m1*m2>= 0){
+          if(((m1>0) && (m2>0)) || ((m1< 0) && (m2< 0))){
             continue;
           }
-          if(fabs(m1)<=eps0){
-            ptd3[0]=p1d3[0];  ptd3[1]=p1d3[1];  ptd3[2]=p1d3[2]; 
+          if(m1==0){
+            ptd3[0]=p1d3[0]; ptd3[1]=p1d3[1]; ptd3[2]=p1d3[2];
             ptuv[0]=puv1[0]; ptuv[1]=puv1[1];
-          }
-          else{
-            ptd3[0]=p2d3[0];  ptd3[1]=p2d3[1];  ptd3[2]=p2d3[2]; 
+          }else{
+            ptd3[0]=p2d3[0]; ptd3[1]=p2d3[1]; ptd3[2]=p2d3[2];
             ptuv[0]=puv2[0]; ptuv[1]=puv2[1];
           }
-        }
-        else{
+        }else{  
           ptd3[0]=1/(m1-m2)*(-m2*p1d3[0]+m1*p2d3[0]);
           ptd3[1]=1/(m1-m2)*(-m2*p1d3[1]+m1*p2d3[1]);
           ptd3[2]=1/(m1-m2)*(-m2*p1d3[2]+m1*p2d3[2]);
-          ptuv[0]=1/(m1-m2)*(-m2*puv1[0]+m1*puv2[1]);
+          ptuv[0]=1/(m1-m2)*(-m2*puv1[0]+m1*puv2[0]);
+          ptuv[1]=1/(m1-m2)*(-m2*puv1[1]+m1*puv2[1]);
         }
         tmpd3[0]=ptd3[0]-pta[0];
-        tmpd3[1]=ptd3[1]-pta[1]; 
+        tmpd3[1]=ptd3[1]-pta[1];
         tmpd3[2]=ptd3[2]-pta[2]; 
-        tmp1d1=dotprod(3,tmpd3,vec1);
-        tmp2d1=fabs(tmp1d1);
-        tmpd3[0]=p1d3[0]-p2d3[0];
-        tmpd3[1]=p1d3[1]-p2d3[1]; 
-        tmpd3[2]=p1d3[2]-p2d3[2]; 
-        tmp3d1=dotprod(3,tmpd3,vec1);
-        tmp4d1=dist(3,p1d3,p2d3);
-        tmp6d1=tmp3d1/tmp4d1;
-        if(tmp1d1>eps0){
-          if(fabs(tmp6d1)>0.75){
-            nptL++;
-            ptL[nptL][0]=ptd3[0]; ptL[nptL][1]=ptd3[1]; ptL[nptL][2]=ptd3[2];
-            ptL[nptL][3]=p2d3[0]; ptL[nptL][4]=p2d3[1]; ptL[nptL][5]=p2d3[2];
-            continue;
-          }
-          if(tmp2d1>epspt){
-            return 1;
-          }
-          if(tmp2d1<2*epspt){
-            continue;
+        crossprod(3,tmpd3,Vec, tmp1d3);
+        if(norm(3,tmp1d3)<Eps1){
+		  tmpd1=zparapt(ptd3)-zparapt(pta)-Eps2;
+          if(tmpd1>0){
+            return tmpd1;
+          }else{
+            out=fmax(out,tmpd1);
           }
         }
       }
     }
   }
-  if(nptL==0){
-    return 0;
-  }
-  for(i=1; i<=nsL; i++){
-    p2d3[0]=sL[i][3]; p2d3[1]=sL[i][4]; p2d3[2]=sL[i][5]; 
-    flg=0;
-    for(j=1; j<=nsL; j++){
-      p1d3[0]=sL[j][0]; p1d3[1]=sL[j][1]; p1d3[2]=sL[j][2]; 
-      if(dist(3,p1d3,p2d3)<eps0){
-        sL[i][6]=j;
-        flg=1;
-        break;
-      }
-    }
-    if((flg==0)&&(sL[i][6]>0)){
-      sL[i][6]=-2;
-    }
-  }
-  for(i=1; i<=nptL; i++){
-    p2d3[0]=ptL[i][3]; p2d3[1]=ptL[i][4]; p2d3[2]=ptL[i][5]; 
-    flg=0;
-    for(j=1; j<=nsL; j++){
-      p1d3[0]=sL[j][0]; p1d3[1]=sL[j][1]; p1d3[2]=sL[j][2];
-      if(dist(3,p1d3,p2d3)<eps0){
-        flg=1;
-        break;
-      }
-    }
-    if(flg==0){
-      ptd3[0]=ptL[i][0]; ptd3[1]=ptL[i][1]; ptd3[2]=ptL[i][2];
-      tmpd1=dist(3,p2d3,ptd3);
-      tmpd3[0]=p2d3[0]-ptd3[0]; 
-      tmpd3[1]=p2d3[1]-ptd3[1];
-      tmpd3[2]=p2d3[2]-ptd3[2]; 
-      tmp1d1=fabs(dotprod(3,tmpd3,vec1));
-      tmp2d1=sqrt(fabs(pow(tmpd1,2.0)-pow(tmp1d1,2.0)));
-      tmp3d1=tmp2d1*dist(3,ptd3,pta)/tmpd1;
-      if(tmp3d1>epspt/5){
-        return 1;
-      }
-      else{
-        continue;
-      }
-    }
-    nxt=sL[j][6];
-    while(nxt>0){
-      ptd3[0]=sL[nxt][3]; ptd3[1]=sL[nxt][4]; ptd3[2]=sL[nxt][5];
-      nxt=sL[nxt][6];
-    }
-    if(nxt==-2){
-      return 1;
-    }
-    else{
-      out=0;
-    }
-  }
   return out;
 }
 
-int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
-  double eps0=pow(10,-5.0), s, p1[2], p2[2], q[2],tmpd1, tmpd2[2], epstmp[2];
-  double fh[DsizeM][2], pta[3], ptap[2], vec[3];
+int nohiddenpara2(short ch,double par[], double fk[][3], int uveq, double figkL[][3]){
+  double eps0=pow(10,-4.0), s, p1[2], p2[2], q[2],tmpd1, tmpd2[2],tmp1,tmp2;
+  double paL[DsizeM], fh[DsizeM][2], pta[3], ptap[2], vec[3],flgd;
   double figL[DsizeM][2], tmpmd2[DsizeM][2], tmpmd3[DsizeM][3];
   double tmp1d2[2], tmp2d2[2], tmp3d2[2], tmpd3[3],tmpd4[4];
   double pd2[2], qd2[2], pd3[3], qd3[3], sp, sq, tp, tq, dtmp, dtmp1, dtmp2;
-  int ncusp=floor(Cusppt[0][0]+0.5), seL[DsizeM], npaL=floor(paL[0]+0.5);
-  int i, n, nfh, nfk, nfigL, nfigkL, nhidden, cspflg, nseL=0, ntmp, flg, tmpi1[DsizeM];
+  int ncusp=floor(Cusppt[0][0]+0.5), seL[DsizeM];
+  int i,j,n, nfh, nfk, cspflg,tmpi1[DsizeM];
+  vec[0]=sin(THETA)*cos(PHI);
+  vec[1]=sin(THETA)*sin(PHI);
+  vec[2]=cos(THETA);
   nfk=length3(fk);
   nfh=projpara(fk, fh);
   pull2(1,fh,p1);
@@ -1007,30 +766,34 @@ int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
      continue;
     }
   }
-  seL[0]=0; nseL=0;
-  for(n=1; n<=npaL-1; n++){
-    s=(paL[n]+paL[n+1])/2.0;
-    ntmp=invparapt(s,fh,fk, pta);
+  paL[0]=0;
+  tmp2=-1;
+  for(j=1;j<=length1(par)-1;j++){
+    tmp1=par[j];
+    if(fabs(tmp1-tmp2)>Eps){
+      add1(paL,tmp1);
+      tmp2=tmp1;
+    }
+  }
+  tmp1=paL[length1(paL)];
+  tmp2=par[length1(par)];
+  if(fabs(tmp1-tmp2)<Eps){
+    paL[length1(paL)]=tmp2;
+  }else{
+    add1(paL,tmp2);
+  }
+  seL[0]=0;
+  for(n=1; n<=length1(paL)-1; n++){
+	s=(paL[n]+paL[n+1])/2.0;
+    invparapt(s,fh,fk, pta);
     parapt(pta,ptap);
-    epstmp[0]=Eps1; epstmp[1]=1;
-    if((n==1)&&(cspflg % 2==0)){
-      pull2(1,fh,p1);
-      dtmp=fmin(Eps1*dist(2,ptap, p1),Eps1); 
-	  epstmp[0]=dtmp; epstmp[1]=1; 
-    }
-    if((n==npaL-1)&&(cspflg % 3==0)){
-      pull2(nfh,fh,p1);
-      dtmp=fmin(Eps1*dist(2, ptap, p1), Eps1); 
-	  epstmp[0]=dtmp; epstmp[1]=1; 
-    }
-    flg=pthiddenQ(Ch,pta, uveq, epstmp);
-    if(flg==0){
-      nseL=appendi1(seL,n);
+    flgd=pthiddenQ(ch,pta, uveq);
+    if(flgd<=0){ //180518
+      appendi1(n,seL);
     }
   }
-  figL[0][0]=0; nfigL=0;
-  figkL[0][0]=0; nfigkL=0;
-  for(i=1; i<=nseL; i++){
+  figL[0][0]=0;figkL[0][0]=0;
+  for(i=1; i<=seL[0]; i++){
     pointoncurve(paL[seL[i]], fh, Eps, tmp1d2);
     pointoncurve(paL[seL[i]+1], fh, Eps, tmp2d2);
     if(i==1){
@@ -1045,12 +808,12 @@ int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
         sq=paL[seL[i]+1];
       }
       else{
-        ntmp=partcrv(sp,sq,fh,tmpmd2);
-        nfigL=appendd2(2,1,ntmp,tmpmd2,figL);
+        partcrv(sp,sq,fh,tmpmd2);
+        appendd2(2,1,length2(tmpmd2),tmpmd2,figL);
         tp=invparapt(sp,fh,fk, tmpd3);
         tq=invparapt(sq,fh,fk,tmpd3);
-        ntmp=partcrv3(tp,tq,fk,tmpmd3);
-        nfigkL=appendd3(2,1,ntmp,tmpmd3,figkL);
+        partcrv3(tp,tq,fk,tmpmd3);
+        appendd3(2,1,length3(tmpmd3),tmpmd3,figkL);
          copyd(0,1,tmp1d2,pd2);
         sp=paL[seL[i]];
         copyd(0,1,tmp2d2,qd2);
@@ -1058,34 +821,34 @@ int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
       }
     }
   }
-  if(nseL>0){
+  if(seL[0]>0){
     if(dist(2,pd2,qd2)>Eps1){
       nearestptpt(pd2,fh,tmpd4);
       dtmp1=tmpd4[2];
       nearestptpt(qd2,fh,tmpd4);
       dtmp2=tmpd4[2];
-      ntmp=partcrv(dtmp1,dtmp2,fh,tmpmd2);
-      nfigL=appendd2(2,1,ntmp,tmpmd2,figL);
+      partcrv(dtmp1,dtmp2,fh,tmpmd2);
+      appendd2(2,1,length2(tmpmd2),tmpmd2,figL);
       tp=invparapt(sp,fh,fk,tmpd3);
       tq=invparapt(sq,fh,fk,tmpd3);
-      ntmp=partcrv3(tp,tq,fk,tmpmd3);
-      nfigkL=appendd3(2,1,ntmp, tmpmd3, figkL);
+      partcrv3(tp,tq,fk,tmpmd3);
+      appendd3(2,1,length3(tmpmd3), tmpmd3, figkL);
     }
     else{
-      nfigL=appendd2(2,1,nfh, fh, figL);
-      nfigkL=appendd3(2,1,nfk,fk,figkL);
+      appendd2(2,1,nfh, fh, figL);
+      appendd3(2,1,nfk,fk,figkL);
     }
   }
-  tmpi1[0]=0; ntmp=0;
-  for(i=1; i<=npaL-1; i++){
+  tmpi1[0]=0;
+  for(i=1; i<=length1(paL)-1; i++){
     if(memberi(i,seL)==0){
-      ntmp=appendi1(tmpi1,i);
+      appendi1(i,tmpi1);
     }
   }
-  copyi(0,ntmp,tmpi1,seL);
-  seL[0]=ntmp; nseL=ntmp;
-  Hiddendata[0][0]=0; nhidden=0;
-  for(i=1; i<=nseL; i++){
+  copyi(0,tmpi1[0],tmpi1,seL);
+  seL[0]=tmpi1[0];
+  Hiddendata[0][0]=0;
+  for(i=1; i<=seL[0]; i++){
     dtmp=paL[seL[i]];
     pointoncurve(dtmp,fh,Eps,tmp1d2);
     dtmp=paL[seL[i]+1];
@@ -1104,8 +867,8 @@ int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
       else{
         tp=invparapt(sp,fh,fk,tmpd3);
         tq=invparapt(sq,fh,fk,tmpd3);
-        ntmp=partcrv3(tp,tq,fk,tmpmd3);
-        nhidden=appendd3(2,1,ntmp,tmpmd3,Hiddendata);
+        partcrv3(tp,tq,fk,tmpmd3);
+        appendd3(2,1,length3(tmpmd3),tmpmd3,Hiddendata);
         copyd(0,1,tmp1d2,pd2);
         sp=paL[seL[i]];
         copyd(0,1,tmp2d2,qd2);
@@ -1113,106 +876,108 @@ int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]){
       }
     }
   }
-  if(nseL>0){
+  if(seL[0]>0){
     if(dist(2,pd2,qd2)>Eps1){
       tp=invparapt(sp,fh,fk,tmpd3);
       tq=invparapt(sq,fh,fk,tmpd3);
-      ntmp=partcrv3(tp,tq,fk,tmpmd3);
-      nhidden=appendd3(2,1,ntmp,tmpmd3,Hiddendata);
+      partcrv3(tp,tq,fk,tmpmd3);
+      appendd3(2,1,length3(tmpmd3),tmpmd3,Hiddendata);
     }
     else{
-      nhidden=appendd3(2,1,nfk,fk,Hiddendata);
+      appendd3(2,1,length3(fk),fk,Hiddendata);
     }
   }
-  return nfigkL;
+  return length3(figkL);
 }
 
-int borderparadata(short ch,double fkL[][3], double fsL[][3]){
+void borderparadata(short ch,double fkL[][3], double fsL[][3]){
   double ekL[DsizeL][3], fall[DsizeL][2], fbdxy[DsizeL][2];
   double tmpd3[DsizeM][3],  tmpd2[DsizeM][2], p[2], p3[3];
   double figkL[DsizeM][3], du, dv;
-  int i, j, nall, nekL, nfsL, nfkL=length3(fkL), nfall, nfbdxy;
-  int ndin, din[DsizeS][2], ndin2,din2[DsizeM][2];
-  int ntmp, ntmpd3,ntmpd2, npar=0, nfigkL;
+  int i, j;
+  int din[DsizeS][2],din2[DsizeM][2];
+  int ntmp, ntmpd3,ntmpd2;
   double par[DsizeM];
   du=(Urng[1]-Urng[0])/Mdv;
   dv=(Vrng[1]-Vrng[0])/Ndv;
   fsL[0][0]=0;
-  ekL[0][0]=0; nekL=0;
+  ekL[0][0]=0;;
   if(DrawE==1){
-    tmpd2[0][0]=0; nall=0;
+    tmpd2[0][0]=0;;
     for(j=0; j<=Ndv; j++){
-      nall=add2(tmpd2, Urng[1], Vrng[0]+j*dv);
+      add2(tmpd2, Urng[1], Vrng[0]+j*dv);
     }
-    nall=surfcurve(ch,2, tmpd2,tmpd3);
-    nekL=appendd3(2,1,nall,tmpd3,ekL);
+    surfcurve(ch,tmpd2,tmpd3);
+    appendd3(2,1,length3(tmpd3),tmpd3,ekL);
   }
   if(DrawN==1){
-    tmpd2[0][0]=0; nall=0;
+    tmpd2[0][0]=0;
     for(j=0; j<=Mdv; j++){
-      nall=add2(tmpd2, Urng[0]+j*du, Vrng[1]);
+      add2(tmpd2, Urng[0]+j*du, Vrng[1]);
     }
-    nall=surfcurve(ch,2, tmpd2,tmpd3);
-    nekL=appendd3(2,1,nall,tmpd3,ekL);
+    surfcurve(ch,tmpd2,tmpd3);
+    appendd3(2,1,length3(tmpd3),tmpd3,ekL);
   }
   if(DrawW==1){
-    tmpd2[0][0]=0; nall=0;
+    tmpd2[0][0]=0;
     for(j=0; j<=Ndv; j++){
-      nall=add2(tmpd2, Urng[0], Vrng[0]+j*dv);
+      add2(tmpd2, Urng[0], Vrng[0]+j*dv);
     }
-    nall=surfcurve(ch,2, tmpd2,tmpd3);
-    nekL=appendd3(2,1,nall,tmpd3,ekL);
+    surfcurve(ch,tmpd2,tmpd3);
+    appendd3(2,1,length3(tmpd3),tmpd3,ekL);
   }
   if(DrawS==1){
-    tmpd2[0][0]=0; nall=0;
+    tmpd2[0][0]=0;
     for(j=0; j<=Mdv; j++){
-      nall=add2(tmpd2, Urng[0]+j*du, Vrng[0]);
+      add2(tmpd2, Urng[0]+j*du, Vrng[0]);
     }
-    nall=surfcurve(ch,2, tmpd2,tmpd3);
-    nekL=appendd3(2,1,nall,tmpd3,ekL);
+    surfcurve(ch,tmpd2,tmpd3);
+    appendd3(2,1,length3(tmpd3),tmpd3,ekL);
   }
-  if(nekL>0){
-    nfkL=appendd3(2,1,nekL,ekL,fkL);
+  if(length3(ekL)>0){
+    appendd3(2,1,length3(ekL),ekL,fkL);
   }
-  nfall=projpara(fkL,fall);
-  nfbdxy=makexybdy(ch,fbdxy);
-  ndin=dataindexd2(2,fall,din);
+  projpara(fkL,fall);
+  makexybdy(ch,fbdxy);
+  dataindexd2(2,fall,din);
   Borderpt[0][0]=0;
   Otherpartition[0]=0;
-  for(i=1; i<=ndin-1; i++){
+  for(i=1; i<=length2i(din)-1; i++){
     Otherpartition[i]=Inf; Otherpartition[0]++;
   }
   Borderhiddendata[0][0]=0;  
-  ndin2=dataindexd3(2,fkL,din2);
-  for(i=1; i<= ndin2; i++){
+  dataindexd3(2,fkL,din2);
+  for(i=1; i<= length2i(din2); i++){
 	tmpd3[0][0]=0;
-    ntmpd3=appendd3(0,din2[i][0],din2[i][1],fkL,tmpd3);
-    ntmpd2=projpara(tmpd3,tmpd2);
-    partitionseg(tmpd2,fall, 1,par);
+    appendd3(0,din2[i][0],din2[i][1],fkL,tmpd3);
+    projpara(tmpd3,tmpd2);
+	partitionseg(tmpd2,fall,0,par);
+/*
 	if(length1(par)>2){
-      ntmp=npar-1;
-      tmpd2[0][0]=0; ntmpd2=0;
+      ntmp=length1(par)-1;
+      tmpd2[0][0]=0;
       for(j=2; j<=ntmp; j++){
         pull2(j,Partitionpt,p);
-        ntmpd2=addptd2(tmpd2, p);
+        addptd2(tmpd2, p);
       }
-      ntmp=appendd2(0, 1,ntmpd2,tmpd2,Borderpt);
+      appendd2(0, 1,length2(tmpd2),tmpd2,Borderpt);
     }
-    tmpd3[0][0]=0;
-    for(j=din2[i][0]; j<=din2[i][1]; j++){
-      pull3(j,fkL,p3);
-      ntmpd3=addptd3(tmpd3,p3);
-    }
-	nfigkL=nohiddenpara2(par,tmpd3,1, figkL);
-	if(nfigkL>0){
-      nfsL=appendd3(2,1,nfigkL,figkL, fsL);
+*/
+	tmpd3[0][0]=0;
+    appendd3(0,din2[i][0],din2[i][1],fkL,tmpd3);
+//    for(j=din2[i][0]; j<=din2[i][1]; j++){
+//      pull3(j,fkL,p3);
+//      addptd3(tmpd3,p3);
+//    }
+	nohiddenpara2(ch,par,tmpd3,1, figkL);
+	if(length3(figkL)>0){
+      appendd3(2,1,length3(figkL),figkL, fsL);
     }
-    ntmp=length3(Hiddendata);
+	ntmp=length3(Hiddendata);
     if(ntmp>0){
       appendd3(2,1,ntmp,Hiddendata,Borderhiddendata);
     }
   }
-  return nfsL;
 }
 
 int dropnumlistcrv3(double qd[][3], double eps1, int out[]){
@@ -1226,48 +991,47 @@ int dropnumlistcrv3(double qd[][3], double eps1, int out[]){
     pd[0][0]=0; npd=0;
     npd=appendd3(0,se,en,qd,pd);
     ptL[0]=0; nptL=0;
-    nptL=appendi1(ptL,1);
+    nptL=appendi1(1,ptL);
     pull3(1,pd,p);
     for(k=2; k<=npd-1; k++){
       pull3(k,pd,tmp2d);
       if(dist(3,p,tmp2d)>eps1){
-        nptL=appendi1(ptL,k);
+        nptL=appendi1(k,ptL);
         pull3(k,pd,p);
       }
     }
     pull3(npd-1,pd,p);
     pull3(npd,pd,tmp2d);
     if(dist(3,p,tmp2d)>eps){  // eps -> eps1 ? 
-      nptL=appendi1(ptL,npd);
+      nptL=appendi1(npd,ptL);
     }
-//    if(nptL==1){ptL[0]=0; nptL=0;}//180318db
     if(nout>0){
-      nout=appendi1(out,Infint);
+      nout=appendi1(Infint,out);
     }
     for(i=1; i<=nptL; i++){
-      nout=appendi1(out,ptL[i]);
+      nout=appendi1(ptL[i],out);
     }
   }
   out[0]=nout;
   return nout;
 }
 
-int sfbdparadata(short ch,double outd3[][3]){
+void sfbdparadata(short ch,double outd3[][3]){
   double pdatad3[DsizeL][3];
   double pts[DsizeL][2],out3md[DsizeL][2],eps;
   double tmpmd[DsizeL][2],tmp1md[DsizeL][2],tmp2md[DsizeL][2],tmp3md[DsizeL][2];
   double tmpd[2],tmp1d[2],tmpd3[3],tmp1d3[3],tmp2d3[3],tmp2md3[DsizeL][3];
   int din[DsizeS][2],nlist[DsizeL];
   int ii,jj,kk,n,nall,flg;
-  envelopedata(ch, tmp3md);
-  out3md[0][0]=0;
+  envelopedata(ch, out3md); //180517from
+  tmp3md[0][0]=0;
   pts[0][0]=0;
-  dataindexd2(2,tmp3md,din);
+  dataindexd2(2,out3md,din);
   for(jj=1;jj<=length2i(din);jj++){
 	tmp1md[0][0]=0;
     tmp2md3[0][0]=0;
     for(kk=din[jj][0];kk<=din[jj][1];kk++){
-      pull2(kk,tmp3md,tmpd);
+      pull2(kk,out3md,tmpd);
       addptd2(tmp1md,tmpd);
       surffun(ch,tmpd[0],tmpd[1],tmpd3);
       addptd3(tmp2md3,tmpd3);
@@ -1283,12 +1047,14 @@ int sfbdparadata(short ch,double outd3[][3]){
         addptd2(tmpmd,tmpd);
       }
       if(length2(tmpmd)>0){
-        appendd2(2,1,length2(tmpmd),tmpmd,out3md);
+        appendd2(2,1,length2(tmpmd),tmpmd,tmp3md);
       }
-    }      
+    }
   }
+  out3md[0][0]=0;
+  appendd2(0,1,length2(tmp3md),tmp3md,out3md);
   tmp3md[0][0]=0;
-  appendd2(2,1,length2(pts),pts,tmp3md);
+  appendd2(0,1,length2(pts),pts,tmp3md);
   pts[0][0]=0;
   for(ii=1;ii<=length2(tmp3md);ii++){
     pull2(ii,tmp3md,tmpd);
@@ -1307,17 +1073,12 @@ int sfbdparadata(short ch,double outd3[][3]){
     }
   }
   for(ii=1;ii<=length2(pts);ii++){
-    tmp3md[0][0]=0;
-    appendd2(2,1,length2(out3md),out3md,tmp3md);
-    out3md[0][0]=0;
     pull2(ii,pts,tmpd);
     surffun(ch,tmpd[0],tmpd[1],tmp1d3);
-    dataindexd2(2,tmp3md,din);
+    dataindexd2(2,out3md,din);
     for(jj=1;jj<=length2i(din);jj++){
       tmp2md[0][0]=0;
-      appendd2(2,din[jj][0],din[jj][1],tmp3md,tmp2md);
-      pull2(jj,tmp3md,tmpd);
-      surffun(ch,tmpd[0],tmpd[1],tmp2d3);
+      appendd2(0,din[jj][0],din[jj][1],out3md,tmp2md);
       eps=1;
       for(kk=1;kk<=length2(tmp2md)-1;kk++){
         pull2(kk,tmp2md,tmpd);
@@ -1331,404 +1092,557 @@ int sfbdparadata(short ch,double outd3[][3]){
       pull2(1,tmp2md,tmpd);
       surffun(ch,tmpd[0],tmpd[1],tmpd3);
       if(dist(3,tmp1d3,tmpd3)<eps+Eps){
-        appendd2(2,ii,ii,pts,out3md);
-        for(kk=1;kk<=length2(tmp2md);kk++){
-          pull2(kk,tmp2md,tmpd);
-          addptd2(out3md,tmpd);
-        }
+        prependd2(ii,ii,pts,tmp2md);
+        replacelistd2(2,jj,out3md,tmp2md);
+        dataindexd2(2,out3md,din);//180517
         continue;
       }
       pull2(length2(tmp2md),tmp2md,tmpd);
       surffun(ch,tmpd[0],tmpd[1],tmpd3);
       if(dist(3,tmp1d3,tmpd3)<eps+Eps){
-        appendd2(2,1,length2(tmp2md),tmp2md,out3md);
-        pull2(ii,pts,tmpd);
-        addptd2(out3md,tmpd);
+        appendd2(0,ii,ii,pts,tmp2md);
+        replacelistd2(2,jj,out3md,tmp2md);
+        dataindexd2(2,out3md,din);//180517
         continue;
       }
     }
-    appendd2(2,1,length2(tmp2md),tmp2md,out3md);
   }
   nall=cuspsplitpara(ch,out3md, pdatad3);
   n=length2(Cuspsplitpt);
   Cusppt[0][0]=0;
   nall=appendd2(0, 1,n,Cuspsplitpt,Cusppt);
-  nall= borderparadata(ch,pdatad3, outd3);
+  borderparadata(ch,pdatad3, outd3);
   push3(Inf,3,0,outd3);
   nall=appendd3(0,1,length3(Borderhiddendata),Borderhiddendata,outd3);
   push3(Inf,3,0,outd3);
-  return nall;
 }
 
-void makevaltablemeet(short ch,double pa[3], double vec[3], double eps0){
-  double u, v, du, dv, pt[3];
-  int i, j;
-  du=(Urng[1]-Urng[0])/Mdv;
-  dv=(Vrng[1]-Vrng[0])/Ndv;
-  for(j=0; j<=Ndv; j++){
-    v=Vrng[0]+j*dv;
-    for(i=0; i<=Mdv; i++){
-      u=Urng[0]+i*du;
-      surffun(ch,u, v, pt);
-      if((fabs(vec[1])>eps0)||(fabs(vec[0])>eps0)){
-        Zval[j][i]=(vec[1]*(pt[0]-pa[0])-vec[0]*(pt[1]-pa[1]));
+void groupnearpt3(double qlist[][3],double plist[][3]){
+// 180428
+  int ii,jj,flg,din[DsizeS][2];
+  double p1[3],p2[3],gl[DsizeS][3];
+  while(length3(qlist)>0){
+    gl[0][0]=0;
+    appendd3(2,1,1,qlist,gl);
+    removed3(1,qlist);
+    for(ii=1;ii<=length3(qlist);ii++){
+      flg=0;
+      pull3(ii,qlist,p1);
+      for(jj=1;jj<=length3(gl);jj++){
+        pull3(jj,gl,p2);
+        if(dist(3,p1,p2)<Eps1){
+          flg=1;
+          break;
+        }
       }
-      else{
-        Zval[j][i]=pt[0]-pa[0];
+      if(flg==1){
+        addptd3(gl,p1);
+        removed3(ii,qlist);
       }
     }
+    appendd3(2,1,length3(gl),gl,plist);
   }
-  for(i=0; i<=Mdv; i++){
-    Xval[i]=Urng[0]+i*du;
-  }
-  for(j=0; j<=Ndv; j++){
-    Yval[j]=Vrng[0]+j*dv;
+}
+
+double funmeet(short ch,double u, double v,double pa[3], double vec[3]){
+  double pt[3], out;
+  surffun(ch,u, v, pt);
+  if((fabs(vec[1])>Eps)||(fabs(vec[0])>Eps)){
+    out=(vec[1]*(pt[0]-pa[0])-vec[0]*(pt[1]-pa[1]));
+  }else{
+    out=pt[0];
   }
+  return out;
 }
 
-int meetpoints(short ch,double pta[3], double ptb[3], double eps, double ptL[][3]){
-  double eps0=pow(10,-4.0), vec[3], out[DsizeM], pL[5][2], vL[5];
-  double qL[DsizeM][2], p1d2[2], p2d2[2], tmpd2[2], puv1d2[2], puv2d2[2];
-  double m1,m2, p1d3[3], p2d3[3], v1, v2, v3, ptd3[3], ptuvd2[2], tmp1d1;
-  double po[3]={0,0,0}, tmpmd1[DsizeM], tmp1d3[3];
-  int nout, nptL, ntmp, nqL, i, j, k;
-  vec[0]=ptb[0]-pta[0];  vec[1]=ptb[1]-pta[1];  vec[2]=ptb[2]-pta[2]; 
-  out[0]=0; nout=0;
-  ptL[0][0]=0; nptL=0;
-  makevaltablemeet(Ch,pta,vec,eps0);
-  for(j=0; j<=Ndv-2; j++){
-    for(i=0; i<=Mdv-2; i++){
-      pL[0][0]=Xval[i]; pL[0][1]=Yval[j]; vL[0]=Zval[j][i];
-      pL[1][0]=Xval[i+1]; pL[1][1]=Yval[j]; vL[1]=Zval[j][i+1];
-      pL[2][0]=Xval[i+1]; pL[2][1]=Yval[j+1]; vL[2]=Zval[j+1][i+1];
-      pL[3][0]=Xval[i]; pL[3][1]=Yval[j+1]; vL[3]=Zval[j+1][i];
-      pL[4][0]=Xval[i]; pL[4][1]=Yval[j]; vL[4]=Zval[j][i];
-	  qL[0][0]=0; nqL=0;
+void meetpoints(short ch,double pta[3], double ptb[3], int uveq,double out[][3]){
+  int i, j, k,flg,din[DsizeS][2];
+  double vec[3],v1,v2,v3;
+  double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
+  double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
+  double pl[5][2], vl[5], ql[11][2],p1[2],p2[2],m1,m2;
+  double ptuv[2], puv1[2], puv2[2],p1d3[3], p2d3[3], ptd3[3];
+  double tmpd1, tmpd2[2],tmpd3[3],tmp1d3[3],tmp1md3[DsizeS][3];
+  vec[0]=ptb[0]-pta[0];  vec[1]=ptb[1]-pta[1];  vec[2]=ptb[2]-pta[2];
+  out[0][0]=0;
+  if(norm(3,vec)<Eps){
+    return;
+  }
+  for(j=0; j<=Ndv-1; j++){
+    vval1=Vrng[0]+j*dv;
+    vval2=Vrng[0]+(j+1)*dv;
+    uval1=Urng[0];
+    eval11=funmeet(ch,uval1,vval1,pta,vec);
+    eval12=funmeet(ch,uval1,vval2,pta,vec);
+    for(i=0; i<=Mdv-1; i++){
+      uval2=Urng[0]+(i+1)*du;
+      eval21=funmeet(ch,uval2,vval1,pta,vec);
+      eval22=funmeet(ch,uval2,vval2,pta,vec);
+	  pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
+      pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
+      pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
+      pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
+	  pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
+      ql[0][0]=0;
       for(k=0; k<=3; k++){
-        pull2(k,pL,p1d2);
-        pull2(k+1,pL,p2d2);
-        m1=vL[k]; m2=vL[k+1];
-		if((m1>=eps0)&&(m2>=eps0)){
+        pull2(k,pl,p1);
+        pull2(k+1,pl,p2);
+        m1=vl[k]; m2=vl[k+1];
+        if(fabs(m1)<Eps){
+          addptd2(ql, p1);
           continue;
         }
-        if((m1<=-eps0)&&(m2<=-eps0)){
+        if(fabs(m2)<Eps){
           continue;
         }
-        if(fabs(m1)<eps0){
-          copyd(0,1,p1d2,tmpd2);
-        }
-        else if(fabs(m2)<eps0){
+        if((m1>0)&&(m2>0)){
           continue;
         }
-        else{
-          tmpd2[0]=1/(m1-m2)*(-m2*p1d2[0]+m1*p2d2[0]);
-          tmpd2[1]=1/(m1-m2)*(-m2*p1d2[1]+m1*p2d2[1]);
+        if((m1< 0)&&(m2< 0)){
+          continue;
         }
-        nqL=addptd2(qL,tmpd2);
+        tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
+        tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
+        addptd2(ql, tmpd2);
       }
-      if(nqL==2){
-        pull2(1, qL, puv1d2);
-        pull2(2, qL, puv2d2);
-        surffun(ch,puv1d2[0], puv1d2[1], p1d3);
-        surffun(ch,puv2d2[0], puv2d2[1], p2d3);
+      uval1=uval2;
+      eval11=eval21;
+      eval12=eval22;
+      if(length2(ql)==2){
+        puv1[0]=ql[1][0]; puv1[1]=ql[1][1]; 
+        puv2[0]=ql[2][0]; puv2[1]=ql[2][1];  
+        surffun(ch,puv1[0],puv1[1],p1d3);
+        surffun(ch,puv2[0],puv2[1],p2d3);
         v1=vec[0]; v2=vec[1]; v3=vec[2];
-        if(fabs(v1)>eps0){
+        if(fabs(v1)>Eps){
           m1=pta[2]+v3/v1*(p1d3[0]-pta[0])-p1d3[2];
           m2=pta[2]+v3/v1*(p2d3[0]-pta[0])-p2d3[2];
-        }
-        else if(fabs(v2)>eps0){
+        }else if(fabs(v2)>Eps){
           m1=pta[2]+v3/v2*(p1d3[1]-pta[1])-p1d3[2];
           m2=pta[2]+v3/v2*(p2d3[1]-pta[1])-p2d3[2];
-        }
-        else{
+        }else{
           m1=pta[1]-p1d3[1];
-          m2=pta[1]-p2d3[1];
+          m2=pta[1]-p2d3[2];
         }
-        if(m1*m2>=-eps0){
-          if(((m1>eps0)&&(m2>eps0)) || ((m1<-eps0)&&(m2<-eps0))){
+        if(m1*m2>= 0){
+          if(((m1>0) && (m2>0)) || ((m1< 0) && (m2< 0))){
             continue;
           }
-          if(fabs(m1)<=eps0){
-            copyd(0,2,p1d3,ptd3);
-            copyd(0,1,puv1d2,ptuvd2);
-          }
-          else{
-            copyd(0,2,p2d3,ptd3);
-            copyd(0,1,puv2d2,ptuvd2);
+          if(m1==0){
+            ptd3[0]=p1d3[0]; ptd3[1]=p1d3[1]; ptd3[2]=p1d3[2];
+            ptuv[0]=puv1[0]; ptuv[1]=puv1[1];
+          }else{
+            ptd3[0]=p2d3[0]; ptd3[1]=p2d3[1]; ptd3[2]=p2d3[2];
+            ptuv[0]=puv2[0]; ptuv[1]=puv2[1];
           }
-        }
-        else{  
+        }else{  
           ptd3[0]=1/(m1-m2)*(-m2*p1d3[0]+m1*p2d3[0]);
           ptd3[1]=1/(m1-m2)*(-m2*p1d3[1]+m1*p2d3[1]);
           ptd3[2]=1/(m1-m2)*(-m2*p1d3[2]+m1*p2d3[2]);
+          ptuv[0]=1/(m1-m2)*(-m2*puv1[0]+m1*puv2[0]);
+          ptuv[1]=1/(m1-m2)*(-m2*puv1[1]+m1*puv2[1]);
         }
-        tmp1d3[0]=ptd3[0]-pta[0]; tmp1d3[1]=ptd3[1]-pta[1];
-        tmp1d3[2]=ptd3[2]-pta[2];
-        tmp1d1=dotprod(3,tmp1d3,vec);
-        tmp1d1=tmp1d1/pow(dist(3,vec,po),2.0);
-        if((tmp1d1>-eps0)&&(tmp1d1<1+eps0)){
-          nout=appendd1(out,tmp1d1);
+        tmpd3[0]=ptd3[0]-pta[0];
+        tmpd3[1]=ptd3[1]-pta[1];
+        tmpd3[2]=ptd3[2]-pta[2]; 
+        crossprod(3,tmpd3,vec, tmp1d3);
+        if(norm(3,tmp1d3)<Eps1){
+          tmpd1=dotprod(3,tmpd3,vec);
+          if((tmpd1>-Eps)&&(norm(3,tmpd3)<norm(3,vec)+Eps)){
+            if(length3(out)==0){
+              addptd3(out,ptd3);
+            }else{
+              pull3(length3(out),out,p1d3);
+              if(dist(3,p1d3,ptd3)>Eps){
+                addptd3(out,ptd3);
+              }
+            }
+          }
         }
       }
     }
   }
-  if(nout>>0){
-    simplesort(out);
-    tmpmd1[1]=out[1]; ntmp=1;
-    for(i=2; i<=nout; i++){
-      tmp1d3[0]=(out[i]-tmpmd1[ntmp])*vec[0];
-      tmp1d3[1]=(out[i]-tmpmd1[ntmp])*vec[1];
-      tmp1d3[2]=(out[i]-tmpmd1[ntmp])*vec[2];
-      tmp1d1=dist(3,tmp1d3,po);
-      if(tmp1d1>eps){
-        ntmp=appendd1(tmpmd1,out[i]);
-      }
-    }
-    nout=ntmp;
-    copyd(0,ntmp,tmpmd1,out);
-    for(i=1; i<=nout; i++){
-      tmp1d3[0]=pta[0]+out[i]*vec[0];
-      tmp1d3[1]=pta[1]+out[i]*vec[1];
-      tmp1d3[2]=pta[2]+out[i]*vec[2];
-      nptL=addptd3(ptL,tmp1d3);
+  tmp1md3[0][0]=0;
+  groupnearpt3(out,tmp1md3);
+  dataindexd3(2,tmp1md3,din);
+  out[0][0]=0;
+  for(i=1;i<=din[0][0];i++){
+    p1d3[0]=0; p1d3[1]=0; p1d3[2]=0;
+    tmpd1=din[i][1]-din[i][0]+1;
+    for(j=din[i][0];j<=din[i][1];j++){
+      p1d3[0]=p1d3[0]+tmp1md3[j][0]/tmpd1;
+      p1d3[1]=p1d3[1]+tmp1md3[j][1]/tmpd1;
+      p1d3[2]=p1d3[2]+tmp1md3[j][2]/tmpd1;
     }
+    addptd3(out,p1d3);
   }
-  return nptL;
 }
 
-int crvsfparadata(double fkL[][3], double fbdkL[][3], int sepflg, double out[][3]){
-  double fbdy[DsizeL][2], fk[DsizeM][3],fkp[DsizeM][3], outh[DsizeL][3];
+void crvsfparadata(short ch, double fkL[][3], double fbdkL[][3], int sepflg, double out[][3]){
+  double fbdy[DsizeL][2], fk[DsizeM][3],outh[DsizeL][3];
+//  double fkp[DsizeM][3], dt;
   double fh[DsizeM][2], parL[DsizeM], tmpmd3[DsizeM][3];
-  double ptL[DsizeM][3], tmpd2[2], tmpd3[3], tmp1d1, tmp2d1, dt;
+  double ptL[DsizeM][3], tmpd2[2], tmpd3[3], tmp1d1, tmp2d1;
   double po[2]={0,0}, epsd2[2]={Eps1,1}, pa[3], pb[3];
-  int nbor=length3(Borderhiddendata), nfkL,nfbdkL, nfbdy, nfk, nall;
-  int ncshidden, npar, ndin, din[DsizeM][2], ndin2,din2[DsizeS][2];
-  int nn, i, j, k, n, nptL, nout;
-  nfkL=length3(fkL);
-  nfbdy=projpara(fbdkL,fbdy);
+  int nbor=length3(Borderhiddendata);
+  int ncshidden,din[DsizeM][2],din2[DsizeS][2];
+  int nn, i, j, k, n;
+  projpara(fbdkL,fbdy);
   out[0][0]=0;
   outh[0][0]=0;
-  ndin=dataindexd3(2,fkL,din);
-  for(nn=1; nn<=ndin; nn++){
-    fkp[0][0]=0;
-    nfk=appendd3(0,din[nn][0],din[nn][1],fkL,fkp);
-    n=fmax(1,floor((fmax(Mdv,Ndv)+1)/(nfk-1)+0.5)); //17.06.16
-    dt=1.0/n;
-	fk[0][0]=0;
-    for(i=1; i<nfk; i++){
-      for(j=0; j<n; j++){
-        tmpd3[0]=fkp[i][0]+(fkp[i+1][0]-fkp[i][0])*j*dt;
-        tmpd3[1]=fkp[i][1]+(fkp[i+1][1]-fkp[i][1])*j*dt;
-        tmpd3[2]=fkp[i][2]+(fkp[i+1][2]-fkp[i][2])*j*dt;
-        addptd3(fk,tmpd3);
-      }
-    }
-    nfk=appendd3(0,nfk,nfk,fkp,fk);
-    nall=projpara(fk,fh);
-    partitionseg(fh, fbdy,1, parL);
-    if(sepflg>=0){
-      for(i=1; i<=nfk-1; i++){
+  dataindexd3(2,fkL,din);
+  for(nn=1; nn<=length2i(din); nn++){
+    fk[0][0]=0;
+    appendd3(0,din[nn][0],din[nn][1],fkL,fk);
+    projpara(fk,fh);
+    partitionseg(fh,fbdy,0,parL);
+    if(sepflg==0){ //180522
+      for(i=1; i<=length3(fk)-1; i++){
         pull3(i,fk,pa);
         pull3(i+1,fk,pb);
-        nptL=meetpoints(Ch,pa,pb, Eps2,ptL);
-        for(j=1; j<=nptL; j++){
-          pull3(j,ptL,tmpd3);
+        meetpoints(ch,pa,pb,1,ptL);
+        for(j=1; j<=length3(ptL); j++){
+	      pull3(j,ptL,tmpd3);
           parapt(tmpd3,tmpd2);
           tmp1d1=paramoncurve(tmpd2,i,fh);
           tmp2d1=Inf;
-          for(k=1; k<=npar; k++){
+          for(k=1; k<=length1(parL); k++){
              tmp2d1=fmin(tmp2d1,fabs(parL[k]-tmp1d1));
           }
           tmpd3[0]=pb[0]-pa[0]; tmpd3[1]=pb[1]-pa[1]; 
           tmpd3[2]=pb[2]-pa[2]; 
           if(tmp2d1*dist(3,tmpd3,po)>Eps1){
-            npar=appendd1(parL,tmp1d1);
+            appendd1(tmp1d1,parL);
           }
         }
         simplesort(parL);
       }
     }
-    nohiddenpara2(parL,fk,1, tmpmd3);
+    nohiddenpara2(ch,parL,fk,1, tmpmd3);
     appendd3(2,1,length3(tmpmd3),tmpmd3,out);
     appendd3(2,1, length3(Hiddendata),Hiddendata,outh);
   }
   push3(Inf,3,0,out); //17.06.17
-  nout=appendd3(0,1,length3(outh),outh,out);
-  push3(Inf,3,0,out); //17.06.17
-  return nout;
+  appendd3(0,1,length3(outh),outh,out);
 }
 
-int wireparadata(double bdyk[][3], double udv[], double vdv[],  int num, 
-             double fsL[][3]){
-  double fbdy[DsizeL][2], udiv[DsizeM], vdiv[DsizeM], fbdxy[DsizeM][2];
-  double fs[DsizeL][3], u, v, d, pdata[DsizeM][3], fkL[DsizeM][3];
-  double fsLh[DsizeLL][3];
-  double tmpmd2[DsizeM][2], fk[DsizeM][3], par[DsizeM], p[2];
-  int nfbdy, nd, i, j, k, n, nu, nv, nfbdxy, nfsL, npdata, nfkL, nfk, npar,ntmp;
-  int nwire, ndin, din[DsizeM][2], nfs, ndin2, din2[DsizeM][2];
-  nfbdy=projpara(bdyk, fbdy);
-  Wirept[0][0]=0;
-  udiv[0]=0; nu=0;
-  vdiv[0]=0; nv=0;
-  n=floor(udv[0]+0.5);
-  if(n==-1){
-    nd=floor(udv[1]+0.5);
-    d=(Urng[1]-Urng[0])/(nd+1);
-    for(i=1; i<=nd; i++){
-      u=Urng[0]+i*d;
-      nu=appendd1(udiv,u);
+void crv3onsfparadata(short ch,double fk[][3], double fbdyd3[][3], double out[][3]){
+  // 180609 debugged(renewed)
+  double fbdy[DsizeL][2],fh[DsizeL][2],fks[DsizeL][3],fhs[DsizeL][2],par[DsizeM];
+  double tmpmd3[DsizeL][3],outh[DsizeL][3];
+  int i,j,din[DsizeS][2];
+  projpara(fbdyd3,fbdy);
+  projpara(fk,fh);
+  out[0][0]=0;
+  outh[0][0]=0;
+  dataindexd2(2,fh,din);
+  for(i=1;i<=din[0][0];i++){
+    fhs[0][0]=0; fks[0][0]=0;
+    appendd2(0,din[i][0],din[i][1],fh,fhs);
+    appendd3(0,din[i][0],din[i][1],fk,fks);
+    tmpmd3[0][0]=0;
+    if(length2(fhs)>1){
+      partitionseg(fhs,fbdy,0, par);
+      nohiddenpara2(ch,par,fks,1,tmpmd3);
+      appendd3(2,1,length3(tmpmd3),tmpmd3,out);
+      appendd3(2,1,length3(Hiddendata),Hiddendata,outh);
+    }else{
+      appendd3(2,1,1,fks,out);
     }
   }
-  else{
-    for(i=1; i<=n; i++){
-      nu=appendd1(udiv,udv[i]);
-    }
+  connectseg3(out, Eps1,tmpmd3);
+  out[0][0]=0;
+  appendd3(0,1,length3(tmpmd3),tmpmd3,out);
+  connectseg3(outh, Eps1,tmpmd3);
+  appendd3(3,1,length3(tmpmd3),tmpmd3,out);
+}
+
+void crv2onsfparadata(short ch,double fh[][2], double fbdyd3[][3], double out[][3]){
+  double fk[DsizeL][3];
+  surfcurve(ch,fh,fk);
+  crv3onsfparadata(ch,fk,fbdyd3,out);
+}
+
+void wireparadata(short ch,double bdyk[][3], double udata[], double vdata[],const char *fname,const char *fnameh){
+  double crv2d[DsizeL][2],out[DsizeL][3],out1[DsizeL][3],out2[DsizeL][3];
+  double du=(Urng[1]-Urng[0])/Mdv;
+  double dv=(Vrng[1]-Vrng[0])/Ndv;
+  int i,j,din[DsizeS][2];
+  short allflg=0;
+  if(strlen(fnameh)==0){
+    allflg=1;
+  }
+  char var[20]="wire3d";
+  char varh[20]="wireh3d";
+  char varnow[20]={'\0'};
+  char varhnow[20]={'\0'};
+  sprintf(var,"%s%d",var,ch);
+  sprintf(varh,"%s%d",varh,ch);
+  char dirfname[256] = {'\0'};
+  char dirfnameh[256] = {'\0'};
+  char varname[256] = {'\0'};
+  char varnameh[256] = {'\0'};
+  sprintf(varname,"%s%d",var,ch);
+  sprintf(varnameh,"%s%d",varh,ch);
+  sprintf(dirfname,"%s%s",Dirname,fname);
+  sprintf(dirfnameh,"%s%s",Dirname,fnameh);
+  if((length1(udata)==0)&&(allflg==0)){
+    out[0][0]=0;
+    output3(1,"w",varname,dirfname,out);
+    output3(1,"w",varnameh,dirfnameh,out);
   }
-  n=floor(vdv[0]+0.5);
-  if(n==-1){
-    nd=floor(vdv[1]+0.5);
-    d=(Vrng[1]-Vrng[0])/(nd+1);
-    for(i=1; i<=nd; i++){
-      v=Vrng[0]+i*d;
-      nv=appendd1(vdiv,v);
+  for(i=1;i<=length1(udata);i++){
+    crv2d[0][0]=0;
+    for(j=0; j<=Ndv; j++){
+      add2(crv2d, udata[i],Vrng[0]+j*dv);
+    }
+    if(allflg==1){out[0][0]=0;}
+	crv2onsfparadata(ch,crv2d,bdyk,out);
+    if(allflg==0){
+      dataindexd3(3,out,din);
+      out1[0][0]=0; out2[0][0]=0;
+      appendd3(0,din[1][0],din[1][1],out,out1);
+      appendd3(0,din[2][0],din[2][1],out,out2);
+      if(i==1){
+        output3(1,"w",varname,dirfname,out1);
+        output3(1,"w",varnameh,dirfnameh,out2);
+      }else{
+        output3(0,"a",varname,dirfname,out1);
+        output3(0,"a",varnameh,dirfnameh,out2);
+      }
+    }else{
+      sprintf(varnow,"%s%s%d",var,"u",i);
+      sprintf(varhnow,"%s%s%d",varh,"u",i);
+      output3h("a",varnow, varhnow,fname,out);
     }
   }
-  else{
-    for(i=1; i<=n; i++){
-      nv=appendd1(vdiv,vdv[i]);
-    }
-  }
-  nfbdxy=makexybdy(Ch,fbdxy);
-  Wirept[0][0]=0;
-  fsL[0][0]=0; nfsL=0;
-  fsLh[0][0]=0;
-  for(i=1; i<=nu; i++){
-    u=udiv[i];
-    d=(Vrng[1]-Vrng[0])/(num*Ndv);
-    tmpmd2[0][0]=0; ntmp=0;
-    for(j=1;j<=(num*Ndv)+1; j++){
-      v=Vrng[0]+(j-1)*d;
-      p[0]=u; p[1]=v;
-      ntmp=addptd2(tmpmd2,p);
-    }
-    nfkL=cuspsplitpara(Ch,tmpmd2, fkL); 
-    ndin=dataindexd3(2,fkL,din);
-    for(j=1; j<=ndin; j++){
-      fk[0][0]=0;
-      nfk=appendd3(0, din[j][0], din[j][1], fkL,fk);
-      ntmp=projpara(fk,tmpmd2);
-	  partitionseg(tmpmd2, fbdy, 1, par);
-	  ntmp=length2(Partitionpt);
-      if(ntmp>2){
-        for(k=2; k<=ntmp; k++){
-           pull2(k, Partitionpt, p);
-           nwire=addptd3(Wirept,p);
-        }
-      }
-      nfs=nohiddenpara2(par,fk,1, fs);
-      n=length3(Hiddendata);
-      ntmp=appendd3(2,1,n, Hiddendata,fsLh);
-      nfsL=appendd3(2,1,nfs,fs,fsL);
-    }
-  }
-  for(i=1; i<= nv; i++){
-    v=vdiv[i];
-    d=(Urng[1]-Urng[0])/(num*Mdv);
-    tmpmd2[0][0]=0; ntmp=0;
-    for(j=1;j<=(num*Mdv)+1; j++){
-      u=Urng[0]+(j-1)*d;
-      p[0]=u; p[1]=v;
-      ntmp=addptd2(tmpmd2,p);
-    }
-    nfkL=cuspsplitpara(Ch,tmpmd2, fkL); 
-    ndin=dataindexd3(2,fkL,din);
-    for(j=1; j<=ndin; j++){
-      fk[0][0]=0;
-      nfk=appendd3(0, din[j][0], din[j][1], fkL,fk);
-      ntmp=projpara(fk,tmpmd2);
-      partitionseg(tmpmd2, fbdy, 1, par);
-	  ntmp=length2(Partitionpt);
-      if(ntmp>2){
-        for(k=2; k<=ntmp; k++){
-           pull2(k, Partitionpt, p);
-           nwire=addptd3(Wirept,p);
-        }
+  for(j=1;j<=length1(vdata);j++){
+    crv2d[0][0]=0;
+    for(i=0; i<=Mdv; i++){
+      add2(crv2d, Urng[0]+i*du,vdata[j]);
+    }
+    if(allflg==1){out[0][0]=0;}
+    crv2onsfparadata(ch,crv2d,bdyk,out);
+    if(allflg==0){
+      dataindexd3(3,out,din);
+      out1[0][0]=0; out2[0][0]=0;
+      appendd3(0,din[1][0],din[1][1],out,out1);
+      appendd3(0,din[2][0],din[2][1],out,out2);
+      if(i<length1(vdata)){
+        output3(0,"a",varname,dirfname,out1);
+        output3(0,"a",varnameh,dirfnameh,out2);
+      }else{
+        output3(0,"a",varname,dirfname,out1);
+        outputend(dirfname);
+        output3(0,"a",varnameh,dirfnameh,out2);
+        outputend(dirfnameh);
       }
-      nohiddenpara2(par,fk,1, fs);
-      appendd3(2,1,length3(Hiddendata), Hiddendata,fsLh);
-      appendd3(2,1,length3(fs),fs,fsL);
+    }else{
+      sprintf(varnow,"%s%s%d",var,"v",j);
+      sprintf(varhnow,"%s%s%d",varh,"v",j);
+	  output3h("a",varnow, varhnow,fname,out);
     }
   }
-  push3(Inf,3,0,fsL);
-  nfsL=appendd3(0,1,length3(fsLh),fsLh,fsL);
-  push3(Inf,3,0,fsL);
-  return nfsL;
+  if(allflg==0){
+    outputend(dirfname);
+    outputend(dirfnameh);
+  }
 }
 
-/*
-double cutfun(int ch, double u, double v){
-  double p[3],val ;
-  surffun(ch,u,v,p);
-  val=p[0]+p[1]-0.5;
-  return val;
+void intersectcrvsf(const char *wa, short chfd,double crv[][3],const char *fname){
+  double pa[3],pb[3],out[DsizeM][3],ptL[DsizeM][3];
+  int i;
+  short chcut;
+  char chc[10];
+  sprintf(chc,"%d",chfd);
+  char var[]="intercrvsf";
+  char dirfname[256] = {'\0'};
+  sprintf(dirfname,"%s%s",Dirname,fname);
+  char varname[256] = {'\0'};
+  sprintf(varname,"%s%s",var,chc);
+  ptL[0][0]=0;
+  for(i=1;i<length3(crv);i++){
+    pull3(i,crv,pa);
+    pull3(i+1,crv,pb);
+    meetpoints(chfd,pa,pb,1,out);
+    appendd3(0,1,length3(out),out,ptL);
+  }
+  output3(1,wa,varname,dirfname,ptL);
 }
-*/
 
-void makevaltablecut(int chcut){
-  double eps0=pow(10,-6.0), u, v, du, dv;
-  int i, j, nx=0, ny=0; ;
-  du=(Urng[1]-Urng[0])/Mdv;
-  dv=(Vrng[1]-Vrng[0])/Ndv;
-  for(j=0; j<=Ndv; j++){
-    v=Vrng[0]+j*dv;
-    for(i=0; i<=Mdv; i++){
-      u=Urng[0]+i*du;
-      Zval[j][i]=cutfun(Ch,chcut, u,v);
+void sfcutparadata(short chfd, short ncut, double fbdy3[][3],const char *fname,const char *fnameh){
+  double Vec[3]={sin(THETA)*cos(PHI),sin(THETA)*sin(PHI),cos(THETA)};
+  double v1=Vec[0],v2=Vec[1],v3=Vec[2];
+  double out[DsizeM][2],out2[DsizeM][2],outd3[DsizeL][3];
+  double ekL[DsizeL][3],out1d3[DsizeL][3],out2d3[DsizeL][3];
+  double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
+  double pl[5][2], vl[5], ql[11][2],diff1,diff2;
+  double tmp1md[DsizeM][2],tmp2md[DsizeM][2],tmpd[2],tmp1d[2],tmp2d[2];
+  double tmp2md3[DsizeM][3],tmpd3[3],tmp1d3[3];
+  double du=(Urng[1]-Urng[0])/Mdv;
+  double dv=(Vrng[1]-Vrng[0])/Ndv;
+  double p1[2],p2[2],m1,m2,q11,q12,q21,q22,a1,a3,b1,b3;
+  int i,j,k,din[DsizeS][2],ns,ne,rmv[DsizeS];
+  short chcut, allflg=0;
+  if(strlen(fnameh)==0){
+    allflg=1;
+  }
+  char var[20]="sfcut3d";
+  char varh[20]="sfcuth3d";
+  sprintf(var,"%s%d",var,chfd);
+  sprintf(varh,"%s%d",varh,chfd);
+  char varnow[20]={'\0'};
+  char varhnow[20]={'\0'};
+  char dirfname[256] = {'\0'};
+  char dirfnameh[256] = {'\0'};
+  char varname[256] = {'\0'};
+  char varnameh[256] = {'\0'};
+  sprintf(varname,"%s%d",var,chfd);
+  sprintf(varnameh,"%s%d",varh,chfd);
+  sprintf(dirfname,"%s%s",Dirname,fname);
+  sprintf(dirfnameh,"%s%s",Dirname,fnameh);
+  for(chcut=1;chcut<=ncut;chcut++){
+    out[0][0]=0; out2[0][0]=0;
+    for(j=0; j<=Ndv-1; j++){
+      vval1=Vrng[0]+j*dv;
+      vval2=Vrng[0]+(j+1)*dv;
+      uval1=Urng[0];
+      eval11=cutfun(chfd, chcut,uval1, vval1);
+      eval12=cutfun(chfd, chcut,uval1, vval2);
+      for(i=0; i<=Mdv-1; i++){
+        uval2=Urng[0]+(i+1)*du;
+        eval21=cutfun(chfd, chcut,uval2, vval1);
+        eval22=cutfun(chfd, chcut,uval2, vval2);
+	    pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
+        pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
+        pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
+        pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
+	    pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
+        ql[0][0]=0;
+        for(k=0; k<=3; k++){
+          pull2(k,pl,p1);
+          pull2(k+1,pl,p2);
+          m1=vl[k]; m2=vl[k+1];
+          if(fabs(m1)<Eps){
+            addptd2(ql, p1);
+            continue;
+          }
+          if(fabs(m2)<Eps){
+            continue;
+          }
+          if((m1>0)&&(m2>0)){
+            continue;
+          }
+          if((m1< 0)&&(m2< 0)){
+            continue;
+          }
+          tmpd[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
+          tmpd[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
+          addptd2(ql, tmpd);
+        }
+        uval1=uval2;
+        eval11=eval21;
+        eval12=eval22;
+        if(length2(ql)==2){
+          q11=ql[1][0];q12=ql[1][1];q21=ql[2][0];q22=ql[2][1];
+          a1=uval1;a3=uval2;b1=vval1;b3=vval2;
+          if(((q11==a1)&&(q21==a1))||((q11==a3)&&(q21==a3))){
+            if(((q21==b1)&&(q22==b1))||((q21==b3)&&(q22==b3))){
+              appendd2(2,1,2,ql,out2);
+            }else{
+              appendd2(2,1,2,ql,out);
+            }
+            continue;
+          }
+          if(((q12==b1)&&(q22==b1))||((q12==b3)&&(q22==b3))){
+            if(((q11==a1)&&(q21==a1))||((q11==a3)&&(q21==a3))){
+              appendd2(2,1,2,ql,out2);
+            }else{
+              appendd2(2,1,2,ql,out);
+            }
+            continue;
+          }
+          appendd2(2,1,2,ql,out);
+        }
+      }
+    }
+    dataindexd2(2,out2,din);
+    while(length2i(din)>0){
+      ns=din[1][0];ne=din[1][1];
+      tmp1md[0][0]=0;
+      appendd2(0,ns,ne,out2,tmp1md);
+      appendd2(2,ns,ne,out2,out);
+      removelistd2(2,1,out2);
+      dataindexd2(2,out2,din);
+      rmv[0]=0;
+      for(j=1;j<=length2i(din);j++){
+        ns=din[j][0];ne=din[j][1];
+        tmp2md[0][0]=0;
+        appendd2(0,ns,ne,out2,tmp2md);
+        pull2(1,tmp2md,tmp2d);pull2(1,tmp1md,tmp1d);
+        diff1=dist(2,tmp1d,tmp2d);
+        pull2(2,tmp2md,tmp2d);pull2(2,tmp1md,tmp1d);
+        diff1=diff1+dist(2,tmp1d,tmp2d);
+        pull2(1,tmp2md,tmp2d);pull2(2,tmp1md,tmp1d);
+        diff2=dist(2,tmp1d,tmp2d);
+        pull2(2,tmp2md,tmp2d);pull2(1,tmp1md,tmp1d);
+        diff2=diff2+dist(2,tmp1d,tmp2d);
+        if((diff1<Eps)||(diff2<Eps)){
+          appendi1(j,rmv);
+          continue;
+        }
+      }
+      for(j=1;j<=rmv[0];j++){
+        removelistd2(2,rmv[j],out2);
+      }
+      dataindexd2(2,out2,din);
+    }
+    connectseg(out, Eps1, out2);
+    dataindexd2(2,out2,din);
+    outd3[0][0]=0;
+    for(i=1;i<=length2i(din);i++){
+      ns=din[i][0];ne=din[i][1];
+      tmp1md[0][0]=0;
+      appendd2(0,ns,ne,out2,tmp1md);
+      tmp2md3[0][0]=0;
+      for(j=1;j<=length2(tmp1md);j++){
+        pull2(j,tmp1md,tmpd);
+        surffun(chfd,tmpd[0],tmpd[1],tmpd3);
+        if(length3(tmp2md3)>0){
+          pull3(length3(tmp2md3),tmp2md3,tmp1d3);
+          if(dist(3,tmpd3,tmp1d3)<Eps1){
+            continue;
+          }
+        }
+        addptd3(tmp2md3,tmpd3);
+      }
+      appendd3(2,1,length3(tmp2md3),tmp2md3,outd3);
+	}
+    crv3onsfparadata(chfd,outd3,fbdy3,ekL);
+    if(allflg==0){
+      dataindexd3(3,ekL,din);
+      out1d3[0][0]=0; out2d3[0][0]=0;
+      appendd3(0,din[1][0],din[1][1],ekL,out1d3);
+      appendd3(0,din[2][0],din[2][1],ekL,out2d3);
+      if(chcut==1){
+        output3(1,"w",varname,dirfname,out1d3);
+        output3(1,"w",varnameh,dirfnameh,out2d3);
+      }else{
+        output3(0,"a",varname,dirfname,out1d3);
+        output3(0,"a",varnameh,dirfnameh,out2d3);
+      }
+    }else{
+      sprintf(varnow,"%s%d",var,chcut);
+      sprintf(varhnow,"%s%d",varh,chcut);
+      output3h("a",varnow, varhnow,fname,ekL);
     }
   }
-  for(i=0; i<=Mdv; i++){
-    Xval[i]=Urng[0]+i*du;
-  }
-  for(j=0; j<=Ndv; j++){
-    Yval[j]=Vrng[0]+j*dv;
+  if(allflg==0){
+    outputend(dirfname);
+    outputend(dirfnameh);
   }
 }
 
-int sfcutdata(int ch, double ekL[][3]){
-  double cutuvL[DsizeL][2], out1[DsizeM][2], out2[DsizeM][3], pt[3];
-  int i, j, nall, ndin, din[DsizeM][2], nekL,nout1,nout2;
-//  makevaltablecut(ch); //180320
-//  nall=implicitplot(cutuvL); //180320
-  ekL[0][0]=0; nekL=0;
-  ndin=dataindexd2(2,cutuvL,din);
-  for(i=1; i<=ndin; i++){
-    out1[0][0]=0;
-    nout1=appendd2(2,din[i][0],din[i][1],cutuvL,out1);
-    out2[0][0]=0;
-    for(j=1; j<=nout1; j++){
-      surffun(ch,out1[j][0], out1[j][1], pt);
-      nout2=addptd3(out2,pt);
-    }
-    nekL=appendd3(2,1,nout2,out2,ekL);
-  }
-  return nekL;
-}
-
-int crv3onsfparadata(double fk[][3], double fbdyd3[][3], double out[][3]){
-  double fbdy[DsizeL][2], fh[DsizeL][2], tmpmd2[DsizeL][2], par[DsizeM];
-  int n, nall, nfbdy, nfbdxy, npar, nout;
-  out[0][0]=0;
-  nfbdy=projpara(fbdyd3,fbdy);
-//  nall=ospline3(2,fk,2,fk1);
-  projpara(fk,fh);
-  partitionseg(fh,fbdy, 1, par);
-  nohiddenpara2(par,fk, 1,out);
-  push3(Inf,3,0,out);
-  nall=appendd3(0,1,length3(Hiddendata),Hiddendata,out);
-  push3(Inf,3,0,out); //17.06.16
-  return nall;
-}
-
 int projcoordcurve(double curve[][3], double out[][3]){
   double   SP=sin(PHI), CP=cos(PHI), ST=sin(THETA), CT=cos(THETA);
   double pt[3], xz, yz, zz;
@@ -2058,41 +1972,80 @@ int skeletondata3(double data[][3], double r00,
   return length3(allres);
 }
 
-/*
-void writesfbd(const char *var, const char *fname, const char *fnameh,double sf[][3]){
-  sfbdparadata(sf);
-  output3(var, "", fname,2,sf);
-  output3(var, "h", fnameh,2,Borderhiddendata);
-}
-
-void writecrv(const char *var, const char *fname, const char *fnameh, 
-                double crv[][3],double sf[][3]){
-  crvsfparadata(crv, sf, 1,Crvdata);
-  output3(var, "", fname,2,Crvdata);
-  output3(var, "h", fnameh,2,Crvsfhiddendata);
-}
-
-void writewire(const char *var, const char *fname, const char *fnameh,
-              double sf[][3], int num, double udv[], double vdv[]){
-  wireparadata(sf, num, udv,vdv, Wiredata);
-  output3(var, "", fname,2,Wiredata);
-  output3(var, "h", fnameh,2,Hiddenwiredata);
-}
-
-void writesfcut(const char *var, const char *fname, const char *fnameh,
-              double crv[][3], double sf[][3]){
-  sfcutdata(crv);
-  crv3onsfparadata(crv, sf, Crvonsfdata);
-  output3(var, "", fname,2,Crvonsfdata);
-  output3(var, "h", fnameh,2,Crv3onsfhiddendata);
-}
-
-void writecrvonsf(const char *var, const char *fname, const char *fnameh, 
-                double crv[][3],double sf[][3]){
-  crv3onsfparadata(crv, sf, Crvonsfdata);
-  output3(var, "", fname,2,Crvonsfdata);
-  output3(var, "h", fnameh,2,Crv3onsfhiddendata);
+void readoutdata3(const char *fname, const char *var, double data[][3]){
+  double x,y,z;
+  float xx;
+  char dstr[256] = {'\0'};
+  char str[10] = {'\0'},tmp[2]={'\0'};
+  int linectr=0, start=0, jj,nn,nctr;
+  FILE *fp;
+  fp=fopen(fname,"r");
+  if(fp==NULL){
+    printf("file not found\n");
+    return;
+  }
+  data[0][0]=0;
+  nn=strlen(var);
+  while( !feof(fp)){
+    fgets(dstr,250,fp);
+    linectr++;
+    if(strncmp(dstr,var,nn)==0){
+      start=linectr;
+    }else{
+      if(start==0){
+        continue;
+      }
+    }
+    if(linectr==start){
+      continue;
+    }
+    if(strncmp(dstr,"sta",3)==0){
+      if(linectr>start+1){
+        add3(data,Inf,2,0);
+      }
+      continue;
+    }
+    if(strncmp(dstr,"end",3)==0){
+      continue;
+    }
+    if(strncmp(dstr,"[",1)!=0){
+      break;
+    }
+    str[0]='\0';
+    nctr=0;
+    for(jj=2;jj<250;jj++){
+      tmp[0]='\0'; sprintf(tmp,"%s%c",tmp,dstr[jj]);
+      if(strncmp(tmp,"/",1)==0){
+        break;
+      }
+      if(strncmp(tmp," ",1)==0){
+        continue;
+      }
+      if(strncmp(tmp,",",1)==0){
+        nctr++;
+        if(nctr==1){
+          x=atof(str); str[0]='\0';
+          continue;
+        }
+        if(nctr==2){
+          y=atof(str);str[0]='\0';
+          continue;
+        }
+      }
+      if(strncmp(tmp,"]",1)==0){
+        z=atof(str);str[0]='\0';
+        add3(data,x,y,z);
+        nctr=0;
+        jj++;
+        tmp[0]='\0'; sprintf(tmp,"%s%c",tmp,dstr[jj]);
+        if(strncmp(tmp,",",1)==0){
+          jj++;
+        }
+        continue;
+      }
+      sprintf(str,"%s%s",str,tmp);
+    }
+  }
+  fclose(fp);
 }
-*/
-
 
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h
index df31a4d812f..0adc7779221 100644
--- a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h
@@ -8,7 +8,7 @@ short Ch=1;
 double Cuspsplitpt[DsizeL][2];
 double Cusppt[DsizeM][2];
 double Otherpartition [DsizeM];
-double Partitionpt[DsizeM][2];
+//double Partitionpt[DsizeM][2];
 double Cuspdata[DsizeM];
 double Wirept[DsizeM][3];
 double Hiddendata[DsizeL][3];
@@ -16,43 +16,41 @@ double Borderpt[DsizeM][2];
 double Borderhiddendata[DsizeL][3];
 double Surf[DsizeL][3];
 double Crvdata[DsizeL][3];
-double Wiredata[DsizeLL][3];
+//double Wiredata[DsizeLL][3];
 
-int output3(const char *var, const char *fname, int level, double out[][3]);
-int output3h(const char *var, const char *varh, const char *fname, double out[][3]);
+int output3(short head,const char *wa,const char *var, const char *fname,double out[][3]);
+int output3h(const char *wa,const char *var, const char *varh, const char *fname, double out[][3]);
+void outputend(const char *fname);
 int writedataC(const char *fname, double out[][3]);
 int readdataC(const char *fname, double data[][3]);
+void readoutdata3(const char *fname, const char *var, double data[][3]);
 void parapt(double pt3[3],double pt2[2]);
+double zparapt(double p[3]);
 int projpara(double p3[][3], double p2[][2]);
 double invparapt(double ph, double fh[][2], double fk[][3],double out[3]);
-int surfcurve(short ch,int level, double crv[][2], double pdata[][3]);
-int xyzax3data(double x0, double x1, double y0, double y1, double z0, double z1, 
-                           double out[][3]);
+void surfcurve(short ch,double crv[][2], double pdata[][3]);
 double evlpfun(short ch, double u, double v);
 int envelopedata(short ch, double out[][2]);
-void evlptablepara(short ch);
-//int implicitplot(double out[][2]);
 int cuspsplitpara(short ch,double pdata[][2], double outkL[][3]);
-double funpthiddenQ(short ch,double u, double v, double pa[3]);
 int makexybdy(short ch,double ehL[][2]);
 void partitionseg(double fig[][2],double gL[][2], int ns, double parL[]);
-void makevaltablepth(short ch,double pta[3]);
-int pthiddenQ(short ch,double pta[3], int uveq, double eps[2]);
-int nohiddenpara2(double paL[], double fk[][3], int uveq, double figkL[][3]);
-int borderparadata(short ch,double fkL[][3], double fsL[][3]);
+double funpthiddenQ(short ch,double u, double v,double pa[3]);
+double pthiddenQ(short ch,double pta[3], int uveq);
+int nohiddenpara2(short ch,double paL[], double fk[][3], int uveq, double figkL[][3]);
+void borderparadata(short ch,double fkL[][3], double fsL[][3]);
 int dropnumlistcrv3(double qd[][3], double eps1, int out[]);
-int sfbdparadata(short ch,double outd3[][3]);
-void makevaltablemeet(short ch,double pa[3], double vec[3], double eps0);
-int meetpoints(short ch,double pta[3], double ptb[3], double eps, double ptL[][3]);
-int crvsfparadata(double fkL[][3], double fbdkL[][3], int sepflg, double out[][3]);
-int sfcutdata(int ch, double ekL[][3]);
-int crv3onsfparadata(double fk[][3], double fbdyd3[][3], double out[][3]);
-int wireparadata(double bdyk[][3],  double udv[], double vdv[], int num, 
-          double fsL[][3]);
+void sfbdparadata(short ch,double outd3[][3]);
+double funmeet(short ch,double u, double v,double pa[3], double vec[3]);
+void meetpoints(short ch,double pta[3], double ptb[3], int uveq,double out[][3]);
+void crvsfparadata(short ch,double fkL[][3], double fbdkL[][3], int sepflg, double out[][3]);
+void crv3onsfparadata(short ch,double fk[][3], double fbdyd3[][3], double out[][3]);
+void crv2onsfparadata(short ch,double fh[][2], double fbdyd3[][3], double out[][3]);
+void intersectcrvsf(const char *wa, short chfd,double crv[][3],const char *fname);
+void wireparadata(short ch,double bdyk[][3], double udata[], double vdata[],const char *fname,const char *fnameh);
+void sfcutparadata(short chfd, short ncut, double fbdy3[][3],const char *fname,const char *fnameh);
 int projcoordcurve(double curve[][3], double out[][3]);
 int kukannozoku(double jokyo[2], double kukanL[][2], double res[][2]);
 int skeletondata3(double data[][3], double r00, double eps1, double eps2,
           double allres[][3]);
 
-double cutfun(short ch,int chcut,double u, double v);
-void makevaltablecut(int chcut);
+double cutfun(short chfd,short chcut,double u, double v);