summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/scripts/ketcindy/ketlibC
diff options
context:
space:
mode:
authorNorbert Preining <preining@logic.at>2018-04-09 16:31:54 +0000
committerNorbert Preining <preining@logic.at>2018-04-09 16:31:54 +0000
commit002ddfda9f7dd612f60f72cc767373c0d0aaf0ec (patch)
treebfe43d0b833c6b6af524d36fa3723d19c2021f75 /Master/texmf-dist/scripts/ketcindy/ketlibC
parent6d863b89d9160515dba7332b23cce2766d7c761f (diff)
ketcindy 20180409
git-svn-id: svn://tug.org/texlive/trunk@47395 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/scripts/ketcindy/ketlibC')
-rw-r--r--Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h2127
-rw-r--r--Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h97
-rw-r--r--Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h2098
-rw-r--r--Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h58
4 files changed, 4380 insertions, 0 deletions
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
new file mode 100644
index 00000000000..27a5434fe41
--- /dev/null
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
@@ -0,0 +1,2127 @@
+//20180308 norm added
+
+int length1(double gdata[]){
+ return floor(gdata[0]+0.5);
+}
+
+int length2(double gdata[][2]){
+ return floor(gdata[0][0]+0.5);
+}
+
+int length3(double gdata[][3]){
+ return floor(gdata[0][0]+0.5);
+}
+
+int length4(double gdata[][4]){
+ return floor(gdata[0][0]+0.5);
+}
+
+int length5(double gdata[][5]){
+ return floor(gdata[0][0]+0.5);
+}
+
+int length6(double gdata[][6]){
+ return floor(gdata[0][0]+0.5);
+}
+
+int length2i(int gdata[][2]){
+ return gdata[0][0];
+}
+
+void copyd(int from, int upto, double q[], double p[]){
+ int i;
+ for(i=from; i<=upto; i++){
+ p[i]=q[i];
+ }
+}
+
+void copyi(int from, int upto, int q[], int p[]){
+ int i;
+ for(i=from; i<=upto; i++){
+ p[i]=q[i];
+ }
+}
+
+int push2(double x, double y, double ptL[][2]){
+ int nall=length2(ptL);
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y;
+ ptL[0][0]=nall;
+ return nall;
+}
+
+int push3(double x, double y, double z, double ptL[][3]){
+ int nall=length3(ptL);
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y; ptL[nall][2]=z;
+ ptL[0][0]=nall;
+ return nall;
+}
+
+int push4(double x, double y, double z, double w, double ptL[][4]){
+ int nall=length4(ptL);
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y;
+ ptL[nall][2]=z; ptL[nall][3]=w;
+ ptL[0][0]=nall;
+ return nall;
+}
+
+int push5(double x, double y, double z, double w, double u, double ptL[][5]){
+ int nall=length5(ptL);
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y;
+ ptL[nall][2]=z; ptL[nall][3]=w; ptL[nall][4]=u;
+ ptL[0][0]=nall;
+ return nall;
+}
+
+int push6(double x, double y, double z, double w, double u, double v,double ptL[][6]){
+ int nall=ptL[0][0];
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y;ptL[nall][2]=z;
+ ptL[nall][3]=w; ptL[nall][4]=u; ptL[nall][5]=v;
+ ptL[0][0]=nall;
+ return nall;
+}
+
+void pull3(int num, double ptL[][3],double pt[3]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2];
+}
+
+void pull2(int num, double ptL[][2],double pt[2]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1];
+}
+
+void pull4(int num, double ptL[][4],double pt[4]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1];
+ pt[2]=ptL[num][2]; pt[3]=ptL[num][3];
+}
+
+void pull5(int num, double ptL[][5],double pt[5]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1];
+ pt[2]=ptL[num][2]; pt[3]=ptL[num][3]; pt[4]=ptL[num][4];
+}
+
+void pull6(int num, double ptL[][6],double pt[6]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2];
+ pt[3]=ptL[num][3]; pt[4]=ptL[num][4]; pt[5]=ptL[num][5];
+}
+
+int add2(double gdata[][2], double x, double y){
+ int nall;
+ nall=push2(x, y, gdata);
+ return nall;
+}
+
+int add3(double gdata[][3], double x, double y, double z){
+ int nall;
+ nall=push3(x, y, z,gdata);
+ return nall;
+}
+
+int addptd3(double gdata[][3], double pt[3]){
+ int nall;
+ nall=push3(pt[0],pt[1],pt[2],gdata);
+ return nall;
+}
+
+int addptd2(double gdata[][2], double pt[2]){
+ int nall;
+ nall=push2(pt[0],pt[1],gdata);
+ return nall;
+}
+
+int addptd4(double gdata[][4], double pt[4]){
+ int nall;
+ nall=push4(pt[0],pt[1],pt[2],pt[3],gdata);
+ return nall;
+}
+
+int addptd5(double gdata[][5], double pt[5]){
+ int nall;
+ nall=push5(pt[0],pt[1],pt[2],pt[3],pt[4],gdata);
+ return nall;
+}
+
+int addptd6(double gdata[][6], double pt[6]){
+ int nall;
+ nall=push6(pt[0],pt[1],pt[2],pt[3],pt[4],pt[5],gdata);
+ return nall;
+}
+
+int appendd2(int level, int from,int upto, double gdata[][2],double mat[][2]){
+ int j, np, nall=length2(mat);
+ np=upto-from+1;
+ if(np<0){np=0;}
+ if(np>0){
+ if(level>0){
+ if(nall>0){
+ nall=add2(mat,Inf, level);
+ }
+ }
+ 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[0][0]=nall+np;
+ return nall+np;
+}
+
+int appendd3(int level, int from,int upto, double gdata[][3],double mat[][3]){
+ int j, np, nall=length3(mat);
+ np=upto-from+1;
+ if(np<0){np=0;}
+ if(np>0){
+ if(level>0){
+ if(nall>0){
+ nall=add3(mat,Inf, level, 0);
+ }
+ }
+ 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 prependd2(int from,int upto, double gdata[][2],double mat[][2]){
+ int j,np, nall=length2(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];
+ }
+ for(j=1; j<=np; j++){
+ mat[j][0]=gdata[from+j-1][0];
+ mat[j][1]=gdata[from+j-1][1];
+ }
+ }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];
+ }
+ 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[0][0]=nall+np;
+ return nall+np;
+}
+
+int appendd1(double numvec[], double num){
+ int n;
+ n=floor(numvec[0]+0.5);
+ numvec[n+1]=num;
+ numvec[0]=n+1;
+ return n+1;
+}
+
+int appendi1(int numvec[],int num){
+ int n;
+ n=numvec[0];
+ numvec[n+1]=num;
+ numvec[0]=n+1;
+ return n+1;
+}
+
+int prependi1(int num, int 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 removei2(int nrmv,int mat[][2]){
+ int j,nall=mat[0][0];
+ 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[0][0]=nall-1;
+ return nall-1;
+ }
+ else{
+ mat[0][0]=0;
+ return 0;
+ }
+}
+
+int dataindexd1(double out[],int din[][2]){
+ int n1=1, j, ctr=0, nall=floor(out[0]+0.5);
+ if(nall==0){
+ din[0][0]=0;
+ return 0;
+ }
+ for(j=1; j<=nall; j++){
+ if(out[j]>Inf-1){
+ 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 dataindexd2(int level,double out[][2],int din[][2]){
+ int n1=1, j, ctr=0, nall=length2(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 dataindexd3(int level, double out[][3],int din[][2]){
+ int n1=1, j, ctr=0, nall=length3(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 dataindexd4(int level, double out[][4],int din[][2]){
+ int n1=1, j, ctr=0, nall=length4(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){
+ din[0][0]=0; din[0][1]=-1;
+ return 0;
+ }
+ for(j=1; j<=nall; j++){
+ if(out[j]==Infint){
+ ctr++;
+ if(j==ninfp+1){
+ din[ctr][0]=0; din[ctr][1]=-1;
+ }
+ else{
+ din[ctr][0]=n1; din[ctr][1]=j-1;
+ }
+ ninfp=j;
+ n1=j+1;
+ }
+ }
+ ctr++;
+ if(out[nall]==Infint){
+ din[ctr][0]=0; din[ctr][1]=0;
+ }
+ else{
+ din[ctr][0]=n1; din[ctr][1]=nall;
+ }
+ din[0][0]=ctr;
+ return ctr;
+}
+
+void dispvec(int n,double dt[]){
+ int i;
+ for(i=0; i<n; i++){
+ printf("%f ",dt[i]);
+ }
+ printf("\n");
+}
+
+void dispveci(int n,int dt[]){
+ int i;
+ for(i=0; i<n; i++){
+ printf("%d ",dt[i]);
+ }
+ printf("\n");
+}
+
+void dispmatd2(int from, int upto, double mat[][2]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %7.5f %7.5f\n",i,mat[i][0],mat[i][1]);
+ }
+}
+
+void dispmatd3(int from, int upto, double mat[][3]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f\n",i,mat[i][0],mat[i][1],mat[i][2]);
+ }
+}
+
+void dispmatd4(int from, int upto, double mat[][4]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f\n",mat[i][2],mat[i][3]);
+ }
+}
+
+void dispmatd5(int from, int upto, double mat[][5]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][2],mat[i][3],mat[i][4]);
+ }
+}
+
+void dispmatd6(int from, int upto, double mat[][6]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f ",i,mat[i][0],mat[i][1],mat[i][2]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][3],mat[i][4],mat[i][5]);
+ }
+}
+
+void dispmatd2all(double mat[][2]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f\n",i,mat[i][0],mat[i][1]);
+ }
+}
+
+void dispmatd3all(double mat[][3]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f\n",i,mat[i][0],mat[i][1],mat[i][2]);
+ }
+}
+
+void dispmatd4all(double mat[][4]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f\n",mat[i][2],mat[i][3]);
+ }
+}
+
+void dispmatd5all(double mat[][5]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][2],mat[i][3],mat[i][4]);
+ }
+}
+
+void dispmatd6all(double mat[][6]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f ",i,mat[i][0],mat[i][1],mat[i][2]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][3],mat[i][4],mat[i][5]);
+ }
+}
+
+void dispmati2(int from, int upto, int mat[][2]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %d %d\n",i,mat[i][0],mat[i][1]);
+ }
+}
+
+void dispmati1(int from, int upto, int mat[]){
+ int i;
+ for(i=from; i<upto+1; i++){
+ printf("%3d %d\n",i,mat[i]);
+ }
+}
+
+void readdata3(const char *fname, double out[][3]){
+ double data[DsizeL][4];
+ int i, j, ndata, nall;
+ 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]);
+ }
+ i++;
+ ndata++;
+ }
+ ndata--;
+ 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]);
+ }
+ return;
+}
+
+void readdata2(const char *fname, double out[][2]){
+ double data[DsizeL][3];
+ int i, j, ndata, nall;
+ 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]);
+ }
+ i++;
+ ndata++;
+ }
+ ndata--;
+ fclose(fp);
+ out[0][0]=0; nall=0;
+ for(i=1; i<=ndata; i++){
+ nall=add2(out, data[i][0], data[i][1]);
+ }
+ return;
+}
+
+int output2(const char *var, const char *fname, int level, double out[][2]){
+ int nall=length2(out), ndin, din[DsizeM][2],i,j;
+ double tmpd[2];
+ FILE *fp;
+ fp=fopen(fname,"w");
+ if (fp == NULL) {
+ printf("cannot open\n");
+ return -1;
+ }
+ ndin=dataindexd2(level, out,din);
+ fprintf(fp,"%s//\n",var);
+ for(i=1; i<=ndin; i++){
+ fprintf(fp,"start//\n");
+ fprintf(fp,"[");
+ for(j=din[i][0]; j<=din[i][1]; j++){
+ pull2(j,out,tmpd);
+ fprintf(fp,"[ %7.5f, %7.5f]",tmpd[0],tmpd[1]);
+ if(j<din[i][1]){
+ fprintf(fp,",");
+ }
+ }
+ fprintf(fp,"]//\n");
+ if(i<ndin){
+ fprintf(fp,"end//\n");
+ }
+ else{
+ fprintf(fp,"end////\n");
+ }
+ }
+ fclose(fp);
+ return 0;
+}
+
+void simplesort(double number[]){
+ int i,j, total=number[0];
+ double tmp;
+ for (i=1; i<=total; i++) {
+ for (j=i+1; j<=total; j++) {
+ if (number[i] > number[j]) {
+ tmp = number[i];
+ number[i] = number[j];
+ number[j] = tmp;
+ }
+ }
+ }
+}
+
+int memberi(int a, int list[]){
+ int i, n=list[0];
+ for(i=1; i<=n; i++){
+ if(a==list[i]){
+ return 1;
+ }
+ }
+ return 0;
+}
+
+double dotprod(int dim, double p[], double q[]){
+ double ans;
+ ans=p[0]*q[0]+p[1]*q[1];
+ if(dim>2){
+ ans=ans+p[2]*q[2];
+ }
+ return ans;
+}
+
+double norm(int dim, double p[]){
+ double ans;
+ ans=dotprod(dim,p,p);
+ return sqrt(ans);
+}
+
+double dist(int dim, double p[], double q[]){
+ double tmp;
+ tmp=pow(q[0]-p[0],2)+pow(q[1]-p[1],2);
+ if(dim>2){
+ tmp=tmp+pow(q[2]-p[2],2);
+ }
+ return sqrt(tmp);
+}
+
+void crossprod(int dim,double a[3],double b[3], double out[3]){
+ if(dim==3){
+ out[0]=a[1]*b[2]-a[2]*b[1];
+ out[1]=a[2]*b[0]-a[0]*b[2];
+ out[2]=a[0]*b[1]-a[1]*b[0];
+ }else{
+ out[0]=a[0]*b[1]-a[1]*b[0];
+ out[1]=Inf;
+ }
+}
+
+void reflectpoint(double p1[2],double regard[][2],double p2[2]){
+ double pa[2],pb[2],u,v,a,b,u2,v2;
+ pull2(1,regard,pa);
+ if(length2(regard)==1){
+ p2[0]=2*pa[0]-p1[0]; p2[1]=2*pa[1]-p1[1];
+ }else{
+ pull2(2,regard,pb);
+ u=pb[0]-pa[0]; u2=pow(u,2.0);
+ v=pb[1]-pa[1]; v2=pow(v,2.0);
+ a=pa[0]; b=pa[1];
+ p2[0]=(u2-v2)/(u2+v2)*p1[0]+2*u*v/(u2+v2)*p1[1];
+ p2[0]=p2[0]-2*v*(u*b-v*a)/(u2+v2);
+ p2[1]=2*u*v/(u2+v2)*p1[0]-(u2-v2)/(u2+v2)*p1[1];
+ p2[1]=p2[1]+2*u*(u*b-v*a)/(u2+v2);
+ }
+}
+
+void pointoncurve(double t, double gdata[][2], double eps, double pt[]){
+ double pa[2],pb[2], s;
+ int n, nall;
+ n=floor(t+eps);
+ s=fmax(t-n,0);
+ nall=length2(gdata);
+ if(n>=nall){
+ pull2(nall,gdata,pt);
+// pt[0]=gdata[nall][0]; pt[1]=gdata[nall][1];
+ }
+ else{
+ pull2(n,gdata,pa);
+ pull2(n+1,gdata,pb);
+ pt[0]=(1-s)*pa[0]+s*pb[0];
+ pt[1]=(1-s)*pa[1]+s*pb[1];
+ }
+}
+
+int partcrv(double a,double b,double pkL[][2], double pL[][2]){
+ int i, is, ie, nall;
+ double eps=pow(10,-3.2);
+ is=ceil(a);
+ ie=floor(b);
+ nall=0;
+ if(a<is-eps){
+ nall++;
+ pL[nall][0]=(is-a)*pkL[is-1][0]+(1-is+a)*pkL[is][0];
+ pL[nall][1]=(is-a)*pkL[is-1][1]+(1-is+a)*pkL[is][1];
+ }
+ for(i=is; i<=ie; i++){
+ nall++;
+ pL[nall][0]=pkL[i][0];
+ pL[nall][1]=pkL[i][1];
+ }
+ if(b>ie+eps){
+ nall++;
+ pL[nall][0]=(1-b+ie)*pkL[ie][0]+(b-ie)*pkL[ie+1][0];
+ pL[nall][1]=(1-b+ie)*pkL[ie][1]+(b-ie)*pkL[ie+1][1];
+ }
+ pL[0][0]=nall;
+ return nall;
+}
+
+int partcrv3(double a,double b,double pkL[][3], double pL[][3]){
+ int i, is, ie, nall;
+ double eps=pow(10,-3.2);
+ is=ceil(a);
+ ie=floor(b);
+ nall=0;
+ if(a<is-eps){
+ nall++;
+ pL[nall][0]=(is-a)*pkL[is-1][0]+(1-is+a)*pkL[is][0];
+ pL[nall][1]=(is-a)*pkL[is-1][1]+(1-is+a)*pkL[is][1];
+ pL[nall][2]=(is-a)*pkL[is-1][2]+(1-is+a)*pkL[is][2];
+ }
+ for(i=is; i<=ie; i++){
+ nall++;
+ pL[nall][0]=pkL[i][0];
+ pL[nall][1]=pkL[i][1];
+ pL[nall][2]=pkL[i][2];
+ }
+ if(b>ie+eps){
+ nall++;
+ pL[nall][0]=(1-b+ie)*pkL[ie][0]+(b-ie)*pkL[ie+1][0];
+ pL[nall][1]=(1-b+ie)*pkL[ie][1]+(b-ie)*pkL[ie+1][1];
+ pL[nall][2]=(1-b+ie)*pkL[ie][2]+(b-ie)*pkL[ie+1][2];
+ }
+ pL[0][0]=nall;
+ return nall;
+}
+
+int connectseg(double pdata[][2], double eps, double out[][2]){
+ int din[DsizeM][2], nd, nq, i, j, flg, st, en, ntmp;
+ int nall=length2(pdata), nout=0;
+ double qd[DsizeM][2], ah[2], ao[2], p[2], q[2], tmp[2];
+ out[0][0]=0;
+ nd=dataindexd2(0,pdata,din);
+ qd[0][0]=0;
+ nq=appendd2(0,din[1][0],din[1][1],pdata,qd);
+ nd=removei2(1,din);
+ while(nd>0){
+ pull2(1,qd,ah);
+ pull2(nq,qd,ao);
+ flg=0;
+ for(j=1; j<=nd; j++){
+ st=din[j][0];
+ en=din[j][1];
+ pull2(st,pdata,p);
+ pull2(en,pdata,q);
+ if(dist(2,p,ao)<eps){
+ nq=appendd2(0,st+1,en,pdata,qd);
+ nd=removei2(j,din);
+ flg=1;
+ break;
+ }
+ if(dist(2,q,ao)<eps){
+ nq=appendd2(0,en-1,st,pdata,qd);
+ nd=removei2(j,din);
+ flg=1;
+ break;
+ }
+ if(dist(2,p,ah)<eps){
+ nq=prependd2(en,st+1,pdata,qd);
+ nd=removei2(j,din);
+ flg=1;
+ break;
+ }
+ if(dist(2,q,ah)<eps){
+ nq=prependd2(st,en-1,pdata,qd);
+ nd=removei2(j,din);
+ flg=1;
+ break;
+ }
+ }
+ if(flg==0){
+ nout=appendd2(2,1,nq,qd,out);
+ qd[0][0]=0;
+ nq=appendd2(0,din[1][0],din[1][1],pdata,qd);
+ nd=removei2(1,din);
+ }
+ }
+ if(nq>0){
+ nout=appendd2(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;
+ double tmp1d2[2], tmp2d1, p[2], q[2];
+ v[0]=b[0]-a[0]; v[1]=b[1]-a[1];
+ sv2=pow(v[0],2.0)+pow(v[1],2.0);
+ if(sv2<pow(10,-6.0)){
+ out[0]=Inf;
+ return;
+ }
+ p[0]=c[0]-a[0]; p[1]=c[1]-a[1];
+ q[0]=d[0]-a[0]; q[1]=d[1]-a[1];
+ epss=fmin(eps2, eps/sqrt(sv2));
+ p1=(p[0]*v[0]+p[1]*v[1])/sv2;
+ p2=(-p[0]*v[1]+p[1]*v[0])/sv2;
+ q1=(q[0]*v[0]+q[1]*v[1])/sv2;
+ q2=(-q[0]*v[1]+q[1]*v[0])/sv2;
+ m1=-epss; M1=1+epss;
+ m2=-epss; M2=epss;
+ if((fmax(p1,q1)<m1) || (fmin(p1,q1)>M1)){
+ out[0]=Inf;
+ return;
+ }
+ if((fmax(p2,q2)<m2) || (fmin(p2,q2)>M2)){
+ out[0]=Inf;
+ return;
+ }
+ if(p2*q2<0){
+ t=p1-(q1-p1)/(q2-p2)*p2;
+ if((t>m1)&& (t<M1)){
+ if((t>-eps0) && (t<1+eps0)){
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ tmp2d1=fmin(fmax(t,0),1);
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=tmp2d1; out[3]=0;
+ }
+ else{
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ tmp2d1=fmin(fmax(t,0),1);
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=tmp2d1; out[3]=1;
+ }
+ return;
+ }
+ if((p1<m1)||(p1>M1)||(p2<m2)||(p2>M2)){
+ if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){
+ out[0]=Inf;
+ return;
+ }
+ t=fmin(fmax(q1,0),1);
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=1;
+ return;
+ }
+ t=fmin(fmax(p1,0),1);
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=1;
+ return;
+ }
+ if((p1>-eps0)&&(p1<1+eps0)&&(p2>-eps0)&&(p2<eps0)){
+ t= p1;
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=0;
+ return;
+ }
+ if((q1>-eps0)&&(q1<1+eps0)&&(q2>-eps0)&&(q2<eps0)){
+ t= q1;
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=0;
+ return;
+ }
+ if((p1<m1)||(p1>M1)||(p2<m2)||(p2>M2)){
+ if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){
+ out[0]=Inf;
+ return;
+ }
+ t=fmin(fmax(q1,0),1);
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=1;
+ return;
+ }
+ if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){
+ t=fmin(fmax(p1,0),1);
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=1;
+ return;
+ }
+ if(fabs(p2)<fabs(q2)){
+ t=fmin(fmax(p1,0),1);
+ }
+ else{
+ t=fmin(fmax(q1,0),1);
+ }
+ tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1];
+ out[0]=tmp1d2[0]; out[1]=tmp1d2[1];
+ out[2]=t; out[3]=1;
+ return;
+}
+
+double paramoncurve(double p[2], int n, double ptL[][2]){
+ double pa[2], pb[2], v[2], w[2], o[2], d2, s;
+ int nptL=length2(ptL);
+ if(n==nptL){
+ return nptL;
+ }
+ else{
+ pull2(n,ptL,pa);
+ pull2(n+1,ptL,pb);
+ v[0]=pb[0]-pa[0]; v[1]=pb[1]-pa[1];
+ w[0]=p[0]-pa[0]; w[1]=p[1]-pa[1];
+ o[0]=0; o[1]=0;
+ d2=pow(dist(2,v,o),2.0);
+ s=dotprod(2,v,w)/d2;
+ s=fmin(fmax(s,0),1);
+ return n+s;
+ }
+}
+
+void nearestptpt(double pa[2], double pL[][2], double pt[4]){
+ double eps0=pow(10,-6.0), p[2], pm[2], im, sm, s;
+ double ans[7], a1,a2, b1,b2, v1, v2, x1, x2, t, tmpd1;
+ int i, n, npL=length2(pL);
+ ans[0]=pa[0]; ans[1]=pa[1];
+ ans[2]=1;
+ p[0]=pL[1][0]; p[1]=pL[1][1];
+ ans[3]=p[0]; ans[4]=p[1];
+ ans[5]=1; ans[6]=dist(2,pa,p);
+ pm[0]=pL[1][0]; pm[1]=pL[1][1];
+ im=1;
+ sm=dist(2,pm,pa);
+ for(i=1; i<=npL-1; i++){
+ a1=pL[i][0]; a2=pL[i][1];
+ b1=pL[i+1][0]; b2=pL[i+1][1];
+ v1=b1-a1; v2=b2-a2;
+ x1=pa[0]; x2=pa[1];
+ tmpd1=pow(v2,2.0)+pow(v1,2.0);
+ if(fabs(tmpd1)<eps0){
+ continue;
+ }
+ t=(-a2*v2-v1*a1+v1*x1+x2*v2)/tmpd1;
+ if(t<-eps0){
+ p[0]=a1; p[1]=a2;
+ }
+ else if(t>1+eps0){
+ p[0]=b1; p[1]=b2;
+ }
+ else{
+ p[0]=a1+t*v1; p[1]=a2+t*v2;
+ }
+ s=dist(2,p,pa);
+ if(s<sm-eps0){
+ tmpd1=paramoncurve(p,i,pL);
+ pm[0]=p[0]; pm[1]=p[1];
+ im=tmpd1; sm=s;
+ }
+ }
+ if(sm<ans[6]){
+ ans[0]=pa[0]; ans[1]=pa[1];
+ ans[2]=n;
+ ans[3]=pm[0]; ans[4]=pm[1];
+ ans[5]=im; ans[6]=sm;
+ }
+ pt[0]=ans[3]; pt[1]=ans[4];
+ pt[2]=ans[5]; pt[3]=ans[6];
+}
+
+int intersectselfPp(double data1[][2],double eps,double eps1s,double kL[][4]){
+ double a[2], b[2], tmpd2[2], tmpd4[4], tmpd5[5], p[2], q[2], t, s;
+ double eps0=pow(10,-5.0);
+ double tmpd1, tmp1d1, tmp2d1;
+ int i, j, k, n, is, ndata1,nkL, ntmp, ntmp1,flg, closed=0;
+ ndata1=length2(data1);
+ p[0]=data1[1][0]; p[1]=data1[1][1];
+ q[0]=data1[ndata1][0]; q[1]=data1[ndata1][1];
+ if(dist(2,p,q)<Eps){closed=1;}
+ kL[0][0]=0; nkL=0;
+ for(i=1; i<=ndata1; i++){
+ ntmp=0; tmpd1=Inf;
+ pull2(i,data1,a);
+ nearestptpt(a, data1, tmpd4);
+ n=floor(tmpd4[2]+0.5);
+ if(n==i){continue;}
+ if(closed==1){
+ if((i==1)&&(n==ndata1)){continue;}
+ if((i==ndata1)&&(n==1)){continue;}
+ }
+ ntmp1=floor(tmpd4[2]+Eps);
+ if((tmpd4[3]<eps1s)&&(ntmp1!=i)){
+ if(ntmp<ntmp1){
+ ntmp=ntmp1;
+ tmpd1=tmpd4[3];
+ nkL=addptd4(kL,tmpd4);
+ }
+ else if(tmpd4[3]<tmpd1){
+ tmpd1=tmpd4[3];
+ kL[0][0]--;
+ nkL=addptd4(kL,tmpd4);
+ }
+ }
+ }
+ return nkL;
+}
+
+/*
+int intersectcrvsPptest(double data1[][2], double data2[][2],
+ double eps1,double eps2, double kL[][4]){
+ double t, s, a[2], b[2], p[2], q[2];
+ double sqeps=pow(10,-10.0), eps0=pow(10,-2.0);
+ double tmp1d2[2],tmp2d2[2],tmpd4[4],tmp1d4[4],tmpd5[5];
+ double tmpd1,tmp1d1,tmp2d1, kL0[DsizeL][5], ds;
+ int i, j, k, n, is, nkL, nall, ndata1,ndata2,ntmp, ntmp1, flg;
+ ndata1=length2(data1);
+ ndata2=length2(data2);
+ if(ndata1==ndata2){
+ tmp1d1=0; tmp2d1=0;
+ for(i=1; i<=ndata1; i++){
+ tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0);
+ tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0);
+ }
+ tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1);
+ if((tmp1d1<Eps)||(tmp2d1<Eps)){
+ kL[0][0]=0; nkL=0;
+ return nkL;
+ }
+ }
+ kL0[0][0]=0; nkL=0;
+ for(i=1; i<=ndata1; i++){
+ ntmp=0; tmpd1=Inf;
+ pull2(i,data1,a);
+ nearestptpt(a, data2, tmpd4);
+ tmp1d1=floor(tmpd4[2]+Eps); //17.06.03
+ a[0]=tmpd4[0]; a[1]=tmpd4[1];
+ nearestptpt(a, data1, tmpd4);
+ tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1];
+ tmpd5[2]=tmpd4[2]; tmpd5[3]=tmp1d1;
+ tmpd5[4]=tmpd4[3];
+ nkL=addptd5(kL0,tmpd5);
+ }
+ kL[0][0]=0; nall=0;
+ p[0]=Inf; p[1]=0;
+ for(i=1; i<=nkL; i++){
+ if(kL0[i][4]<eps2){
+ q[0]=kL0[i][0]; q[1]=kL0[i][1];
+ if(dist(2,p,q)>eps2){
+ nall=push4(kL0[i][0],kL0[i][1],kL0[i][2],kL0[i][3],kL);
+ p[0]=kL0[i][0]; p[1]=kL0[i][1];
+ ds=kL0[i][4];
+ continue;
+ }
+ if(kL0[i][4]<ds){
+ kL[0][0]--;
+ nall=push4(kL0[i][0],kL0[i][1],kL0[i][2],kL0[i][3],kL);
+ p[0]=kL0[i][0]; p[1]=kL0[i][1];
+ ds=kL0[i][4];
+ }
+ }
+ }
+ return nall;
+}
+*/
+
+int intersectcrvsPp(double data1[][2], double data2[][2],
+ double eps,double eps1s, double kL[][4]){
+ double t, s, a[2], b[2], p[2], q[2];
+ double sqeps=pow(10,-10.0), eps0=pow(10,-2.0);
+ double tmp1d2[2],tmp2d2[2],tmpd4[4],tmp1d4[4],tmpd5[5];
+ double tmpd1,tmp1d1,tmp2d1;
+ int i, j, k, n, is, ndata1, ndata2, nkL,nkL1,nkL2,ntmp, ntmp1, flg;
+ ndata1=length2(data1);
+ ndata2=length2(data2);
+ if(ndata1==ndata2){
+ tmp1d1=0; tmp2d1=0;
+ for(i=1; i<=ndata1; i++){
+ tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0);
+ tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0);
+ }
+ tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1);
+ if((tmp1d1<Eps)||(tmp2d1<Eps)){
+ kL[0][0]=0; nkL=0;
+ return nkL;
+ }
+ }
+ kL[0][0]=0; nkL=0;
+ for(i=1; i<=ndata1;i++){
+ ntmp=0; tmpd1=Inf;
+ pull2(i,data1,a);
+ nearestptpt(a, data2, tmpd4);
+ tmp1d1=floor(tmpd4[2]+Eps); //17.06.03
+ a[0]=tmpd4[0]; a[1]=tmpd4[1];
+ nearestptpt(a, data1, tmpd4);
+ if(tmpd4[3]<eps1s){
+ flg=0;
+ p[0]=tmpd4[0]; p[1]=tmpd4[1];
+ for(j=1; j<=nkL; j++){
+ q[0]=kL[j][0]; q[1]=kL[j][1];
+ if(dist(2,p,q)<eps0){
+ flg=1;
+ break;
+ }
+ }
+ if(flg==1){continue;}
+ ntmp1=floor(tmpd4[2]+Eps);
+ if(ntmp<ntmp1){
+ ntmp=ntmp1;
+ tmpd1=tmpd4[3];
+ tmp1d4[0]=tmpd4[0]; tmp1d4[1]=tmpd4[1];
+ tmp1d4[2]=tmpd4[2]; tmp1d4[3]=tmp1d1;
+ nkL=addptd4(kL,tmp1d4);
+ }
+ else if(tmpd4[3]<tmpd1){
+ tmpd1=tmpd4[3];
+ kL[0][0]--;
+ tmp1d4[0]=tmpd4[0]; tmp1d4[1]=tmpd4[1];
+ tmp1d4[2]=tmpd4[2]; tmp1d4[3]=tmp1d1;
+ nkL=addptd4(kL,tmp1d4);
+ }
+ }
+ }
+ return nkL;
+}
+
+int intersectcrvsPpold(double data1[][2], double data2[][2],
+ double eps,double eps1s, double kL[][4]){
+// 17.05.22
+ double kL1[DsizeL][5], kL2[DsizeL][5];
+ double t, s, a[2], b[2], p[2], q[2];
+ double sqeps=pow(10,-10.0), eps0=pow(10,-6.0);
+ double tmp1d2[2],tmp2d2[2],tmpd4[4],tmpd5[5];
+ double tmpd1,tmp1d1,tmp2d1;
+ int i, j, k, n, is, ndata1, ndata2, nkL,nkL1,nkL2,ntmp,flg;
+ ndata1=length2(data1);
+ ndata2=length2(data2);
+ if(ndata1==ndata2){
+ tmp1d1=0; tmp2d1=0;
+ for(i=1; i<=ndata1; i++){
+ tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0);
+ tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0);
+ tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0);
+ }
+ tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1);
+ if((tmp1d1<eps0)||(tmp2d1<eps0)){
+ kL[0][0]=0; nkL=0;
+ return nkL;
+ }
+ }
+ kL1[0][0]=0; nkL1=0;
+ kL2[0][0]=0; nkL2=0;
+ for(i=1; i<=ndata1-1; i++){
+ pull2(i,data1,a);
+ pull2(i+1,data1,b);
+ for(j=1; j<=ndata2-1; j++){
+ pull2(j, data2, p);
+ pull2(j+1, data2, q);
+ koutenseg(a,b,p,q,eps,eps1s,tmpd4);
+ if(tmpd4[0]<Inf-1){
+ ntmp=floor(tmpd4[3]+Eps);
+ if(ntmp==0){
+ tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1];
+ tmpd5[2]=tmpd4[2]; tmpd5[3]=i; tmpd5[4]=j;
+ nkL1=addptd5(kL1,tmpd5);
+ }
+ else{
+ tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1];
+ tmpd5[2]=tmpd4[2]; tmpd5[3]=i; tmpd5[4]=j;
+ nkL2=addptd5(kL2,tmpd5);
+ }
+ }
+ }
+ }
+ kL[0][0]=0; nkL=0;
+ if(nkL1>0){
+ pull5(1,kL1,tmpd5);
+ p[0]=tmpd5[0]; p[1]=tmpd5[1];
+ t=tmpd5[2];
+ i=floor(tmpd5[3]+Eps);
+ j=floor(tmpd5[4]+Eps);
+ tmpd4[0]=p[0]; tmpd4[1]=p[1];
+ tmpd4[2]=i+t; tmpd4[3]=j;
+ nkL=addptd4(kL,tmpd4);
+ }
+ for(n=2; n<=nkL1; n++){
+ pull5(n,kL1,tmpd5);
+ p[0]=tmpd5[0]; p[1]=tmpd5[1];
+ flg=0;
+ for(k=1; k<=nkL; k++){
+ pull4(k,kL,tmpd4);
+ tmp1d2[0]=tmpd4[0]; tmp1d2[1]=tmpd4[1];
+ if(pow(dist(2,p,tmp1d2),2.0)<sqeps){
+ flg=1;
+ break;
+ }
+ }
+ if(flg==0){
+ tmpd4[0]=p[0]; tmpd4[1]=p[1];
+ t=tmpd5[2];
+ i=floor(tmpd5[3]+Eps);
+ j=floor(tmpd5[4]+Eps);
+ tmpd4[2]=t+i; tmpd4[3]=j;
+ nkL=addptd4(kL,tmpd4);
+ }
+ }
+ for(n=1; n<=nkL2; n++){
+ pull5(n,kL2,tmpd5);
+ p[0]=tmpd5[0]; p[1]=tmpd5[1];
+ flg=0;
+ for(k=1; k<=nkL; k++){
+ pull4(k,kL,tmpd4);
+ tmp1d2[0]=tmpd4[0]; tmp1d2[1]=tmpd4[1];
+ if(pow(dist(2,p,tmp1d2),2.0)<sqeps){
+ flg=1;
+ break;
+ }
+ }
+ if(flg==0){
+ tmpd4[0]=p[0]; tmpd4[1]=p[1];
+ t=tmpd5[2];
+ i=floor(tmpd5[3]+Eps);
+ j=floor(tmpd5[4]+Eps);
+ tmpd4[2]=i+t; tmpd4[3]=j;
+ nkL=addptd4(kL,tmpd4);
+ }
+ }
+ return nkL;
+}
+
+int dropnumlistcrv(double qd[][2], double eps1, int out[]){
+ int i,j,k, nout,nall=length2(qd), nd, se, en, npd, nptL;
+ int din[DsizeM][2],ptL[DsizeL];
+ double eps=pow(10.0,-6.0), pd[DsizeL][2], p[2], tmp2d[2];
+ out[0]=0; nout=0;
+ nd=dataindexd2(2,qd,din);
+ for(j=1; j<=nd; j++){
+ se=din[j][0]; en=din[j][1];
+ pd[0][0]=0; npd=0;
+ npd=appendd2(0,se,en,qd,pd);
+ ptL[0]=0; nptL=0;
+ nptL=appendi1(ptL,1);
+ 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);
+ 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);
+ }
+ if(nptL==1){ptL[0]=0; nptL=0;}
+ if(nout>0){
+ nout=appendi1(out,Infint);
+ }
+ for(i=1; i<=nptL; i++){
+ nout=appendi1(out,ptL[i]);
+ }
+ }
+ out[0]=nout;
+ return nout;
+}
+
+int increasept2(double ptL[][2], int nn, double out[][2]){
+ int ii, jj, npt, nall;
+ double tmpd2[2];
+ npt=length2(ptL);
+ out[0][0]=0;
+ for(ii=1; ii<=npt-1; ii++){
+ for(jj=0; jj<nn; jj++){
+ tmpd2[0]=ptL[ii][0]+1.0*jj/nn*(ptL[ii+1][0]-ptL[ii][0]);
+ tmpd2[1]=ptL[ii][1]+1.0*jj/nn*(ptL[ii+1][1]-ptL[ii][1]);
+ addptd2(out,tmpd2);
+ }
+ }
+ nall=appendd2(0,npt,npt,ptL,out);
+ return nall;
+}
+
+int increasept3(double ptL[][3], int nn, double out[][3]){
+ int ii, jj, npt, nall;
+ double tmpd3[3];
+ npt=length3(ptL);
+ out[0][0]=0;
+ for(ii=1; ii<=npt-1; ii++){
+ for(jj=0; jj<nn; jj++){
+ tmpd3[0]=ptL[ii][0]+1.0*jj/nn*(ptL[ii+1][0]-ptL[ii][0]);
+ tmpd3[1]=ptL[ii][1]+1.0*jj/nn*(ptL[ii+1][1]-ptL[ii][1]);
+ tmpd3[2]=ptL[ii][2]+1.0*jj/nn*(ptL[ii+1][2]-ptL[ii][2]);
+ addptd3(out,tmpd3);
+ }
+ }
+ nall=appendd3(0,npt,npt,ptL,out);
+ return nall;
+}
+
+void bezierpt(double t, double p0[], double p3[], double p1[], double p2[],
+ double pout[]){
+ double p4[2], p5[2], p6[2], p7[2], p8[2];
+ int i, flg3=1;
+ if(p2[0]>Inf-1){
+ for(i=0; i<2; i++){
+ p2[i]=p3[i];
+ }
+ flg3=0;
+ }
+ for(i=0; i<2; i++){
+ p4[i]=(1-t)*p0[i]+t*p1[i];
+ }
+ for(i=0; i<2; i++){
+ p5[i]=(1-t)*p1[i]+t*p2[i];
+ }
+ for(i=0; i<2; i++){
+ p6[i]=(1-t)*p2[i]+t*p3[i];
+ }
+ for(i=0; i<2; i++){
+ p7[i]=(1-t)*p4[i]+t*p5[i];
+ }
+ if(flg3==0){
+ for(i=0; i<2; i++){
+ pout[i]=p7[i];
+ }
+ }
+ else{
+ for(i=0; i<2; i++){
+ p8[i]=(1-t)*p5[i]+t*p6[i];
+ }
+ for(i=0; i<2; i++){
+ pout[i]=(1-t)*p7[i]+t*p8[i];
+ }
+ }
+}
+
+int bezier(double ptL[][2], double ctrL[][4], int num, double out[][2]){
+ double p0[2], p3[2], p1[2], p2[2], pt[2], dt=1.0/num;
+ int i, j, n, nall;
+ out[0][0]=0; nall=0;
+ nall=appendd2(0,1,1,ptL,out);
+ n=length2(ptL);
+ for(i=1; i<=n-1; i++){
+ pull2(i,ptL, p0);
+ pull2(i+1,ptL, p3);
+ p1[0]=ctrL[i][0]; p1[1]=ctrL[i][1];
+ p2[0]=ctrL[i][2]; p2[1]=ctrL[i][3];
+ for(j=1; j<=num; j++){
+ bezierpt(j*dt,p0,p3,p1,p2,pt);
+ nall=push2(pt[0],pt[1], out);
+ }
+ }
+ return nall;
+}
+
+int ospline(int level, double ptL[][2], int num, double out[][2]){
+ double ctrlist[DsizeM][4], eps0=pow(10,-6.0), tmp, cc;
+ double p[2], q[2], p0[2], p3[2], p1[2], p2[2], pQ[2], pR[2];
+ double ptlist[DsizeM][2], tmpmd2[DsizeM][2];
+ int nd, i, j, npt, nall, nctr, flg, ndin, closed, din[DsizeM][2], ntmp;
+ out[0][0]=0;
+ if(num==1){
+ npt=length2(ptL);
+ for(i=0; i<=npt; i++){
+ out[i][0]=ptL[i][0]; out[i][1]=ptL[i][1];
+ }
+ nall=length2(out);
+ return nall;
+ }
+ ndin=dataindexd2(level,ptL,din);
+ for(nd=1; nd<=ndin; nd++){
+ ptlist[0][0]=0;
+ npt=appendd2(0,din[nd][0],din[nd][1],ptL,ptlist);
+ if(npt==2){
+ appendd2(2, 1,npt,ptlist,out);
+ continue;
+ }
+ closed=0;
+ pull2(1,ptlist,p);
+ pull2(npt,ptlist,q);
+ if(dist(2,p,q)<eps0){
+ closed=1;
+ }
+ ctrlist[0][0]=0; nctr=0;
+ for(i=1; i<=npt-1; i++){
+ flg=0;
+ pull2(i, ptlist,p1);
+ pull2(i+1, ptlist,p2);
+ if((i==1)||(i==npt-1)){
+ flg=1;
+ if(closed==0){
+ pQ[0]=(p1[0]+2*p2[0])/3.0;
+ pQ[1]=(p1[1]+2*p2[1])/3.0;
+ pR[0]=(2*p1[0]+p2[0])/3.0;
+ pR[1]=(2*p1[1]+p2[1])/3.0;
+ nctr++;
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1];
+ ctrlist[nctr][2]=pR[0]; ctrlist[nctr][3]=pR[1];
+ flg=2;
+ }
+ else{
+ if(i==1){
+ pull2(npt-1, ptlist,p0);
+ pull2(i+2, ptlist,p3);
+ }
+ else{
+ pull2(i-1, ptlist, p0);
+ pull2(2, ptlist, p3);
+ }
+ }
+ }
+ if(flg<=1){
+ if(flg==0){
+ pull2(i-1, ptlist, p0);
+ pull2(i+2, ptlist, p3);
+ }
+ p[0]=p2[0]-p0[0]; p[1]=p2[1]-p0[1];
+ q[0]=p3[0]-p1[0]; q[1]=p3[1]-p1[1];
+ tmp=1+sqrt((1+dotprod(2,p,q)/dist(2,p0,p2)/dist(2,p1,p3))/2);
+ cc=4*dist(2,p1,p2)/3/(dist(2,p0,p2)+dist(2,p1,p3))/tmp;
+ pQ[0]=p1[0]+cc*(p2[0]-p0[0]); pQ[1]=p1[1]+cc*(p2[1]-p0[1]);
+ pR[0]=p2[0]+cc*(p1[0]-p3[0]); pR[1]=p2[1]+cc*(p1[1]-p3[1]);
+ nctr++;
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1];
+ ctrlist[nctr][2]=pR[0]; ctrlist[nctr][3]=pR[1];
+ }
+ }
+ if(closed==0){
+ p1[0]=ctrlist[2][0]; p1[1]=ctrlist[2][1];
+ pull2(2, ptlist, p2);
+ pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]);
+ ctrlist[1][0]=pQ[0]; ctrlist[1][1]=pQ[1];
+ ctrlist[1][2]=Inf; ctrlist[1][3]=0;
+ p1[0]=ctrlist[nctr-1][2]; p1[1]=ctrlist[nctr-1][3];
+ pull2(nctr,ptlist,p2);
+ pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]);
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1];
+ ctrlist[nctr][2]=Inf; ctrlist[1][3]=0;
+ }
+ ctrlist[0][0]=nctr;
+ if(length2(ptlist)>2){
+ ntmp=bezier(ptlist,ctrlist,num, tmpmd2);
+ nall=appendd2(2,1,ntmp,tmpmd2,out);
+ }
+ }
+ return nall;
+}
+
+void bezierpt3(double t, double p0[], double p3[], double p1[], double p2[],
+ double pout[]){
+ double p4[3], p5[3], p6[3], p7[3], p8[3];
+ int i, flg3=1;
+ if(p2[0]>Inf-1){
+ for(i=0; i<3; i++){
+ p2[i]=p3[i];
+ }
+ flg3=0;
+ }
+ for(i=0; i<3; i++){
+ p4[i]=(1-t)*p0[i]+t*p1[i];
+ }
+ for(i=0; i<3; i++){
+ p5[i]=(1-t)*p1[i]+t*p2[i];
+ }
+ for(i=0; i<3; i++){
+ p6[i]=(1-t)*p2[i]+t*p3[i];
+ }
+ for(i=0; i<3; i++){
+ p7[i]=(1-t)*p4[i]+t*p5[i];
+ }
+ if(flg3==0){
+ for(i=0; i<3; i++){
+ pout[i]=p7[i];
+ }
+ }
+ else{
+ for(i=0; i<3; i++){
+ p8[i]=(1-t)*p5[i]+t*p6[i];
+ }
+ for(i=0; i<3; i++){
+ pout[i]=(1-t)*p7[i]+t*p8[i];
+ }
+ }
+}
+
+int bezier3(double ptL[][3], double ctrL[][6], int num, double out[][3]){
+ double p0[3], p3[3], p1[3], p2[3], pt[3], dt=1.0/num;
+ int i, j, n, nall;
+ out[0][0]=0; nall=0;
+ nall=appendd3(0,1,1,ptL,out);
+ n=length3(ptL);
+ for(i=1; i<=n-1; i++){
+ pull3(i,ptL, p0);
+ pull3(i+1,ptL, p3);
+ p1[0]=ctrL[i][0]; p1[1]=ctrL[i][1]; p1[2]=ctrL[i][2];
+ p2[0]=ctrL[i][3]; p2[1]=ctrL[i][4]; p2[2]=ctrL[i][5];
+ for(j=1; j<=num; j++){
+ bezierpt3(j*dt,p0,p3,p1,p2,pt);
+ nall=push3(pt[0],pt[1],pt[2], out);
+ }
+ }
+ return nall;
+}
+
+int ospline3(int level, double ptL[][3], int num, double out[][3]){
+ double ctrlist[DsizeM][6], eps0=pow(10,-6.0), tmp, cc;
+ double p[3], q[3], p0[3], p3[3], p1[3], p2[3], pQ[3], pR[3];
+ double ptlist[DsizeM][3], tmpmd2[DsizeM][3];
+ int nd, i, j, n, npt, nall, nctr, closed, flg, ndin, din[DsizeM][2], ntmp;
+ out[0][0]=0; nall=0;
+ if(num==1){
+ npt=length3(ptL);
+ for(i=0; i<=npt; i++){
+ out[i][0]=ptL[i][0]; out[i][1]=ptL[i][1]; out[i][2]=ptL[i][2];
+ }
+ nall=length3(out);
+ return nall;
+ }
+ ndin=dataindexd3(2,ptL,din);
+ for(nd=1; nd<=ndin; nd++){
+ ptlist[0][0]=0;
+ npt=din[nd][1]-din[nd][0];
+ if(npt<=3){
+ n=appendd3(0, din[nd][0], din[nd][0], ptL, ptlist);
+ for(i=1; i<=npt-1; i++){
+ pull3(din[nd][0]+i,ptL,p);
+ pull3(din[nd][0]+i+1,ptL,q);
+ for(j=1; j<=num; j++){
+ p1[0]=p[0]+1.0*j/num*(q[0]-p[0]);
+ p1[1]=p[1]+1.0*j/num*(q[1]-p[1]);
+ p1[2]=p[2]+1.0*j/num*(q[2]-p[2]);
+ n=addptd3(ptlist,p1);
+ }
+ }
+ nall=appendd3(2, 1,n,ptlist,out);
+ continue;
+ }
+ appendd3(0,din[nd][0],din[nd][1],ptL,ptlist);
+ closed=0;
+ pull3(1,ptlist,p);
+ pull3(npt,ptlist,q);
+ if(dist(3,p,q)<eps0){closed=1;}
+ ctrlist[0][0]=0; nctr=0;
+ for(i=1; i<=npt-1; i++){
+ flg=0;
+ pull3(i, ptlist,p1);
+ pull3(i+1, ptlist,p2);
+ if((i==1)||(i==npt-1)){
+ flg=1;
+ if(closed==0){
+ pQ[0]=(p1[0]+2*p2[0])/3.0;
+ pQ[1]=(p1[1]+2*p2[1])/3.0;
+ pQ[2]=(p1[2]+2*p2[2])/3.0;
+ pR[0]=(2*p1[0]+p2[0])/3.0;
+ pR[1]=(2*p1[1]+p2[1])/3.0;
+ pR[2]=(2*p1[2]+p2[2])/3.0;
+ nctr++;
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2];
+ ctrlist[nctr][3]=pR[0]; ctrlist[nctr][4]=pR[1]; ctrlist[nctr][5]=pR[2];
+ flg=2;
+ }
+ else{
+ if(i==1){
+ pull3(npt-1, ptlist,p0);
+ pull3(i+2, ptlist,p3);
+ }
+ else{
+ pull3(i-1, ptlist, p0);
+ pull3(2, ptlist, p3);
+ }
+ }
+ }
+ if(flg<=1){
+ if(flg==0){
+ pull3(i-1, ptlist, p0);
+ pull3(i+2, ptlist, p3);
+ }
+ p[0]=p2[0]-p0[0]; p[1]=p2[1]-p0[1]; p[2]=p2[2]-p0[2];
+ q[0]=p3[0]-p1[0]; q[1]=p3[1]-p1[1]; q[2]=p3[2]-p1[2];
+ tmp=1+sqrt((1+dotprod(3,p,q)/dist(3,p0,p2)/dist(3,p1,p3))/2);
+ cc=4*dist(3,p1,p2)/3/(dist(3,p0,p2)+dist(3,p1,p3))/tmp;
+ pQ[0]=p1[0]+cc*(p2[0]-p0[0]); pQ[1]=p1[1]+cc*(p2[1]-p0[1]);
+ pQ[2]=p1[2]+cc*(p2[2]-p0[2]);
+ pR[0]=p2[0]+cc*(p1[0]-p3[0]); pR[1]=p2[1]+cc*(p1[1]-p3[1]);
+ pR[2]=p2[2]+cc*(p1[2]-p3[2]);
+ nctr++;
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2];
+ ctrlist[nctr][3]=pR[0]; ctrlist[nctr][4]=pR[1]; ctrlist[nctr][5]=pR[2];
+ }
+ }
+ if(closed==0){
+ p1[0]=ctrlist[2][0]; p1[1]=ctrlist[2][1]; p1[2]=ctrlist[2][2];
+ pull3(2, ptlist, p2);
+ pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]);
+ pQ[2]=p2[2]+(p2[2]-p1[2]);
+ ctrlist[1][0]=pQ[0]; ctrlist[1][1]=pQ[1]; ctrlist[1][2]=pQ[2];
+ ctrlist[1][3]=Inf; ctrlist[1][4]=0; ctrlist[1][5]=0;
+ p1[0]=ctrlist[nctr-1][3]; p1[1]=ctrlist[nctr-1][4]; p1[2]=ctrlist[nctr-1][5];
+ pull3(nctr,ptlist,p2);
+ pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]);
+ pQ[2]=p2[2]+(p2[2]-p1[2]);
+ ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2];
+ ctrlist[nctr][3]=Inf; ctrlist[1][4]=0; ctrlist[1][5]=0;
+ }
+ ctrlist[0][0]=nctr;
+ if(length3(ptlist)>2){
+ ntmp=bezier3(ptlist,ctrlist,num, tmpmd2);
+ nall=appendd3(2,1,ntmp,tmpmd2,out);
+ }
+ }
+ return nall;
+}
+
+void segplot(double x1,double y1,double x2,double y2,double seg[][2]){
+ seg[0][0]=0;
+ push2(x1,y1,seg);
+ push2(x2,y2,seg);
+}
+
+void addsegplot(double x1,double y1,double seg[][2]){
+ push2(x1,y1,seg);
+}
+
+void circledata(double cx,double cy,double r, int num,double out[][2]){
+ double t,dt,pt[2];
+ int i;
+ dt=2*M_PI/num;
+ out[0][0]=0;
+ for(i=0; i<=num; i++){
+ push2(cx+r*cos(i*dt),cy+r*sin(i*dt),out);
+ }
+}
+
+void intersectline(double p1[2],double v1[2], double p2[2],double v2[2],double out[4]){
+ double tmp,d,dt,ds,t,s,tmp1[2],tmp2[2],tmp3[2],pt[2];
+ tmp=dotprod(2,v1,v2);
+ tmp1[0]=tmp; tmp1[1]=pow(norm(2,v1),2.0);
+ tmp2[0]=-pow(norm(2,v2),2.0);tmp2[1]=-tmp;
+ pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1];
+ tmp3[0]=dotprod(2,pt,v2);tmp3[1]=dotprod(2,pt,v1);
+ crossprod(2,tmp1,tmp2,pt); d=pt[0];
+ crossprod(2,v1,v2,pt); tmp=fabs(pt[0]);
+ if(tmp>Eps){
+ crossprod(2,tmp3,tmp2,pt); dt=pt[0];
+ crossprod(2,tmp1,tmp3,pt); ds=pt[0];
+ t=dt/d;
+ s=ds/d;
+ pt[0]=p1[0]+v1[0]*t; pt[1]=p1[1]+v1[1]*t;
+ out[0]=pt[0];out[1]=pt[1];out[2]=t; out[3]=s;
+ }else{
+ pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1];
+ crossprod(2,pt,v1,tmp1);
+ tmp=tmp1[0]/norm(2,v1);
+ out[0]=fabs(tmp); out[1]=Inf;
+ }
+}
+
+void intersectseg(double seg1[][2],double seg2[][2], double out[5]){
+ double p1[2],q1[2],p2[2],q2[2],pt[2],pt1[2],nv[2],pts[DsizeS][4];
+ double q1p1[2],q2p2[2],p1q2[2],q1p2[2],p1p2[2],q1q2[2];
+ double tmp[2],tmp1[2],tmp2[2],tmp3[2],tmp4[2];
+ double tmpd1,x,y,t,s,t1,s1,distance,tmpd[4],tmp3d[7];
+ int ctr=0, i;
+ pull2(1,seg1,p1); pull2(2,seg1,q1);
+ pull2(1,seg2,p2); pull2(2,seg2,q2);
+ q1p1[0]=q1[0]-p1[0];q1p1[1]=q1[1]-p1[1];
+ q2p2[0]=q2[0]-p2[0];q2p2[1]=q2[1]-p2[1];
+ p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1];
+ q1p2[0]=q1[0]-p2[0];q1p2[1]=q1[1]-p2[1];
+ p1p2[0]=p1[0]-p2[0];p1p2[1]=p1[1]-p2[1];
+ p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1];
+ q1q2[0]=q1[0]-q2[0];q1q2[1]=q1[1]-q2[1];
+ if((norm(2,q1p1)<Eps)||(norm(2,q2p2)<Eps)){
+ out[0]=-1; out[1]=Inf;
+ }else{
+ intersectline(p1,q1p1,p2,q2p2,tmpd);
+ if(tmpd[1]<Inf-Eps){
+ pt[0]=tmpd[0];pt[1]=tmpd[1];t=tmpd[2];s=tmpd[3];
+ if((t*(t-1)<Eps)&&(s*(s-1)<Eps)){
+ out[0]=0;out[1]=pt[0];out[2]=pt[1];out[3]=t;out[4]=s;
+ }else{
+ if(t<0){t=0;}; if(t>1){t=1;};
+ if(s<0){s=0;}; if(s>1){s=1;};
+ tmp3d[ctr]=norm(2,p1p2);ctr++;
+ tmp3d[ctr]=norm(2,p1q2);ctr++;
+ tmp3d[ctr]=norm(2,q1p2);ctr++;
+ tmp3d[ctr]=norm(2,q1q2);ctr++;
+ tmp1[0]=q2p2[1];tmp1[1]=-q2p2[0];
+ intersectline(p1,tmp1,p2,q2p2,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,p1);ctr++;
+ }
+ }
+ intersectline(q1,tmp1,p2,q2p2,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,q1);ctr++;
+ }
+ }
+ tmp1[0]=q1p1[1];tmp1[1]=-q1p1[0];
+ intersectline(p2,tmp1,p1,q1p1,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,p2);ctr++;
+ }
+ }
+ intersectline(q2,tmp1,p1,q1p1,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,q2);ctr++;
+ }
+ }
+ tmpd1=tmp3d[0];
+ for(i=1;i<ctr;i++){
+ if(tmp3d[i]<tmpd1){tmpd1=tmp3d[i];};
+ }
+ out[0]=tmpd1;out[1]=pt[0];out[2]=pt[1];
+ out[3]=t;out[4]=s;
+ }
+ }else{
+ distance=tmpd[0];
+ nv[0]=q1p1[1]/norm(2,q1p1);nv[1]=-q1p1[0]/norm(2,q1p1);
+ pts[0][0]=0;
+ intersectline(p1,nv,p2,q2p2,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1];
+ push4(x,y,0,s1,pts);
+ }
+ intersectline(q1,nv,p2,q2p2,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1];
+ push4(x,y,1,s1,pts);
+ }
+ intersectline(p2,nv,p1,q1p1,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1];
+ push4(x,y,s1,0,pts);
+ }
+ intersectline(q2,nv,p1,q1p1,tmpd);
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1];
+ push4(x,y,s1,2,pts);
+ }
+ if(length4(pts)==0){
+ tmpd1=norm(2,p1p2);
+ if(norm(2,p1q2)<tmpd1){tmpd1=norm(2,p1q2);}
+ if(norm(2,q1p2)<tmpd1){tmpd1=norm(2,q1p2);}
+ if(norm(2,q1q2)<tmpd1){tmpd1=norm(2,q1q2);}
+ out[0]=tmpd1;out[1]=Inf;
+ }else{
+ if(distance>Eps1){
+ out[0]=distance;out[1]=Inf;
+ }else{
+ x=0; y=0;
+ for(i=1; i<=length4(pts); i++){
+ x=x+pts[i][0];
+ y=y+pts[i][1];
+ }
+ tmp3[0]=x/length4(pts); tmp3[1]=y/length4(pts);
+ nearestptpt(tmp3,seg1,tmpd);
+ t=tmpd[2];
+ nearestptpt(tmp3,seg2,tmpd);
+ s=tmpd[2];
+ out[0]=distance; out[1]=tmp3[0]; out[2]=tmp3[1];
+ out[3]=t;out[4]=s;
+ }
+ }
+ }
+ }
+}
+
+int osplineseg(double ptlist[5][2],int num,double out[][2]){
+ double p0[2],p1[2],p2[2],p3[2],p2p0[2],p3p1[2],p2p1[2];
+ double tmp,cc,pq[2],pr[2],knots[3][2],ctrlist[2][4];
+ pull2(1,ptlist,p0); pull2(2,ptlist,p1);
+ pull2(3,ptlist,p2); pull2(4,ptlist,p3);
+ p2p0[0]=p2[0]-p0[0];p2p0[1]=p2[1]-p0[1];
+ p3p1[0]=p3[0]-p1[0];p3p1[1]=p3[1]-p1[1];
+ p2p1[0]=p2[0]-p1[0];p2p1[1]=p2[1]-p1[1];
+ tmp=norm(2,p2p0)*norm(2,p3p1);
+ tmp=1+sqrt((1+dotprod(2,p2p0,p3p1)/tmp)/2);
+ cc=4*norm(2,p2p1)/3/(norm(2,p2p0)+norm(2,p3p1))/tmp;
+ pq[0]=p1[0]+cc*p2p0[0]; pq[1]=p1[1]+cc*p2p0[1];
+ pr[0]=p2[0]-cc*p3p1[0];pr[1]=p2[1]-cc*p3p1[1];
+ ctrlist[0][0]=0;
+ push4(pq[0],pq[1],pr[0],pr[1],ctrlist);
+ knots[0][0]=0;
+ addptd2(knots,p1);addptd2(knots,p2);
+ bezier(knots,ctrlist,num,out);
+ return length2(out);
+}
+
+void intersectpartseg(double crv1[][2],double crv2[][2],int ii, int jj,int num,double out[6]){
+ double distance=10*Eps2,seg1[3][2],seg2[3][2],tmpd[2],tmp1d[2],tmp2d[2];
+ double tmp,snang,dst,tmpd5[5],os1[DsizeS][2],os2[DsizeS][2];
+ double p0[2],p1[2],p2[2],p3[2],ptlist[5][2],regard[3][2];
+ double tmp1md[20][2],tmpd3[3],tmp3md3[20][3],tmp2md3[20][3];
+ double tmpd4[4],Eps00=pow(10,-8.0);
+ double sx,sy;
+ int kk,ll,nn;
+ out[0]=Inf;
+ for(kk=1;kk<6;kk++){out[kk]=0;}
+ seg1[0][0]=0;seg2[0][0]=0;
+ appendd2(1,ii,ii+1,crv1,seg1);
+ appendd2(1,jj,jj+1,crv2,seg2);
+ pull2(2,seg1,tmp1d);
+ pull2(1,seg1,tmpd);
+ tmp1d[0]=tmp1d[0]-tmpd[0]; tmp1d[1]=tmp1d[1]-tmpd[1];
+ pull2(2,seg2,tmp2d);
+ pull2(1,seg2,tmpd);
+ tmp2d[0]=tmp2d[0]-tmpd[0]; tmp2d[1]=tmp2d[1]-tmpd[1];
+ crossprod(2,tmp1d,tmp2d,tmpd);
+ snang=fabs(tmpd[0])/(norm(2,tmp1d)*norm(2,tmp2d));
+ intersectseg(seg1,seg2,tmpd5);
+ dst=tmpd5[0];
+ if(dst<Eps){
+ if(dst>-Eps){
+ out[0]=tmpd5[1];out[1]=tmpd5[2];out[2]=ii+tmpd5[3];
+ out[3]=jj+tmpd5[4];out[4]=dst;out[5]=snang;
+ }
+ }else{
+ if(dst<Eps2){
+ os1[0][0]=0;os2[0][0]=0;
+ pull2(1,seg1,tmp1d);
+ pull2(2,seg1,tmp2d);
+ tmpd[0]=tmp2d[0]-tmp1d[0];tmpd[1]=tmp2d[1]-tmp1d[1];
+ if((length2(crv1)==2)||(norm(2,tmpd)>distance-Eps)){
+ appendd2(1,1,2,seg1,os1);
+ }else{
+ pull2(1,seg1,p1);pull2(2,seg1,p2);
+ if(ii==1){
+ pull2(3,crv1,p3);
+ tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p3,regard,p0);
+ }else{
+ if(ii==length2(crv1)-1){
+ pull2(ii-1,crv1,p0);
+ tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p0,regard,p3);
+ }else{
+ pull2(ii-1,crv1,p0);pull2(ii+2,crv1,p3);
+ }
+ }
+ ptlist[0][0]=0;
+ addptd2(ptlist,p0);addptd2(ptlist,p1);
+ addptd2(ptlist,p2);addptd2(ptlist,p3);
+ osplineseg(ptlist,num,os1);
+ }
+ pull2(1,seg2,tmp1d);
+ pull2(2,seg2,tmp2d);
+ if((length2(crv2)==2)||(dist(2,tmp2d,tmp1d)>distance-Eps)){
+ appendd2(1,1,2,seg2,os2);
+ }else{
+ pull2(1,seg2,p1);pull2(2,seg2,p2);
+ if(jj==1){
+ pull2(3,crv2,p3);
+ tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p3,regard,p0);
+ }else{
+ if(jj==length2(crv2)-1){
+ pull2(jj-1,crv2,p0);
+ tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p0,regard,p3);
+ }else{
+ pull2(jj-1,crv2,p0);pull2(jj+2,crv2,p3);
+ }
+ }
+ ptlist[0][0]=0;
+ addptd2(ptlist,p0);addptd2(ptlist,p1);
+ addptd2(ptlist,p2);addptd2(ptlist,p3);
+ osplineseg(ptlist,num,os2);
+ }
+ tmp2md3[0][0]=0;
+ for(kk=1;kk<length2(os1);kk++){
+ for(ll=1;ll<length2(os2);ll++){
+ seg1[0][0]=0; seg2[0][0]=0;
+ pull2(kk,os1,tmpd);addptd2(seg1,tmpd);
+ pull2(kk+1,os1,tmpd);addptd2(seg1,tmpd);
+ pull2(ll,os2,tmpd);addptd2(seg2,tmpd);
+ pull2(ll+1,os2,tmpd);addptd2(seg2,tmpd);
+ intersectseg(seg1,seg2,tmpd5);
+ if((tmpd5[0]<Eps1)&&(tmpd5[1]<Inf-Eps)){
+ if(tmpd5[0]<dst+Eps00){
+ dst=tmpd5[0];
+ tmp3md3[0][0]=0;
+ for(nn=1;nn<=length3(tmp2md3);nn++){
+ pull3(nn,tmp2md3,tmpd3);
+ if(tmpd3[0]<dst+Eps00){
+ addptd3(tmp3md3,tmpd3);
+ }
+ }
+ tmp2md3[0][0]=0;
+ for(nn=1;nn<=length3(tmp3md3);nn++){
+ pull3(nn,tmp3md3,tmpd3);
+ addptd3(tmp2md3,tmpd3);
+ }
+ push3(dst,tmpd5[1],tmpd5[2],tmp2md3);
+ }
+ }
+ }
+ }
+ if(length3(tmp2md3)>0){
+ tmp1md[0][0]=0;
+ sx=0; sy=0;
+ for(nn=1;nn<=length3(tmp2md3);nn++){
+ pull3(nn,tmp2md3,tmpd3);
+ push2(tmpd3[1],tmpd3[2],tmp1md);
+ sx=sx+tmpd3[1];
+ sy=sy+tmpd3[2];
+ }
+ sx=sx/length3(tmp2md3);
+ sy=sy/length3(tmp2md3);
+ p0[0]=sx; p0[1]=sy;
+ pull2(ii,crv1,p1); pull2(ii+1,crv1,p2);
+ tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1];
+ tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0];
+ intersectline(p0,tmpd,p1,tmp1d,tmpd4);
+ tmp=tmpd4[3];
+ if(tmp<0){tmp=0;}
+ if(tmp>1){tmp=1;}
+ out[0]=sx;out[1]=sy;out[2]=ii+tmp;
+ pull2(jj,crv2,p1); pull2(jj+1,crv2,p2);
+ tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1];
+ tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0];
+ intersectline(p0,tmpd,p1,tmp1d,tmpd4);
+ tmp=tmpd4[3];
+ if(tmp<0){tmp=0;}
+ if(tmp>1){tmp=1;}
+ out[3]=jj+tmp;out[4]=dst;out[5]=snang;
+ }
+ }
+ }
+}
+
+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];
+ double s1,e1,s2,e2,dst,tmp1,tmp2,tmp1md6[40][6];
+ rL[0][0]=0;
+ if(length6(ptdL)==0){
+ return;
+ }
+ for(ii=2;ii<=length6(ptdL);ii++){
+ pull6(ii,ptdL,tmp1d6);addptd6(rL,tmp1d6);
+ }
+ tmp1md6[0][0]=0;
+ pull6(1,ptdL,tmp1d6);
+ addptd6(tmp1md6,tmp1d6);
+ kk=tmp1d6[2];
+ if(tmp1d6[2]<kk+Eps00){
+ s1=kk-1-Eps00;e1=s1+2+2*Eps00;
+ }else{
+ s1=kk-Eps00; e1=s1+1+2*Eps00;
+ }
+ kk=tmp1d6[3];
+ if(tmp1d6[3]<kk+Eps00){
+ s2=kk-1-Eps00;e2=s2+2+2*Eps00;
+ }else{
+ s2=kk-Eps00; e2=s2+1+2*Eps00;
+ }
+ ctr=1;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++;
+ }
+ }
+ dst=100;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmp1d6);
+ if(tmp1d6[4]<dst){
+ dst=tmp1d6[4];
+ }
+ }
+ ptdL[0][0]=0;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmp1d6);
+ if(tmp1d6[4]<dst+Eps00){
+ addptd6(ptdL,tmp1d6);
+ }
+ }
+ for(ii=0;ii<=length6(rL);ii++){
+ diffL[ii]=0;
+ }
+ for(ii=1;ii<=ctr;ii++){
+ kk=numL[ii];
+ diffL[kk]=1;
+ }
+ nr=length6(rL);
+ rL[0][0]=0;
+ for(ii=1;ii<=nr;ii++){
+ if(diffL[ii]==0){
+ pull6(ii,rL,tmp1d6);addptd6(rL,tmp1d6);
+ }
+ }
+}
+
+void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,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 tmpd[2],tmp1d[2],tmp2d[2],tmpd6[6],tmp1d6[6];
+ double sx,sy,dst,ans[4];
+ int ii,jj,kk,nall,self,js,je;
+ crv1[0][0]=0;
+ pull2(1,crv1s,tmpd);addptd2(crv1,tmpd);
+ for(ii=2;ii<=length2(crv1s);ii++){
+ pull2(length2(crv1),crv1,tmp1d);
+ pull2(ii,crv1s,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>Eps){
+ addptd2(crv1,tmp2d);
+ }
+ }
+ crv2[0][0]=0;
+ pull2(1,crv2s,tmpd);
+ addptd2(crv2,tmpd);
+ for(ii=2;ii<=length2(crv2s);ii++){
+ pull2(length2(crv2),crv2,tmp1d);
+ pull2(ii,crv2s,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>Eps){
+ addptd2(crv2,tmp2d);
+ }
+ }
+ if(length2(crv1)!=length2(crv2)){
+ self=0;
+ }else{
+ self=1;
+ for(ii=1;ii<=length2(crv1);ii++){
+ pull2(ii,crv1,tmp1d);pull2(ii,crv2,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>0){
+ self=0;
+ break;
+ }
+ }
+ }
+ tmp1md6[0][0]=0;
+ for(ii=1;ii<=length2(crv1)-1;ii++){
+ if(self==0){
+ js=1; je=length2(crv2)-1;
+ }else{
+ js=ii+2; je=length2(crv2)-1;
+ }
+ for(jj=js;jj<=je;jj++){
+ intersectpartseg(crv1,crv2,ii,jj,num,tmpd6);
+ if(tmpd6[0]<Inf-Eps){
+ tmpd[0]=tmpd6[0];tmpd[1]=tmpd6[1];
+ if(length6(tmp1md6)==0){
+ addptd6(tmp1md6,tmpd6);
+ }else{
+ pull6(length6(tmp1md6),tmp1md6,tmp1d6);
+ tmp1d[0]=tmp1d6[0];tmp1d[1]=tmp1d6[1];
+ if(dist(2,tmp1d,tmpd)>Eps1){
+ addptd6(tmp1md6,tmpd6);
+ }
+ }
+ if(self==1){
+ push6(tmpd6[0],tmpd6[1],tmpd6[3],tmpd6[2],tmpd6[4],tmpd6[5],tmp1md6);
+ }
+ }
+ }
+ }
+ out[0][0]=0;
+ 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];
+ }
+ }
+ 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];
+ }
+ 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);
+ }
+ }
+}
+
+void intersectcurves(double crv1[][2],double crv2[][2],int num,double out[][2]){
+ double out6[DsizeS][6],tmpd[2],tmpd6[6];
+ int jj;
+ intersectcurvesPp(crv1,crv2,num,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
new file mode 100644
index 00000000000..7e83a2b9e83
--- /dev/null
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommonhead.h
@@ -0,0 +1,97 @@
+int length1(double gdata[]);
+int length2(double gdata[][2]);
+int length3(double gdata[][3]);
+int length4(double gdata[][4]);
+int length5(double gdata[][5]);
+int length6(double gdata[][6]);
+int length2i(int gdata[][2]);
+void copyd(int from, int upto, double q[], double p[]);
+void copyi(int from, int upto, int q[], int p[]);
+int addptd2(double gdata[][2], double pt[2]);
+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 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]);
+int push3(double x, double y, double z, double ptL[][3]);
+int push4(double x, double y, double z, double w, double ptL[][4]);
+int push5(double x, double y, double z, double w, double u, double ptL[][5]);
+int push6(double x, double y, double z, double w, double u, double v,double ptL[][6]);
+void pull2(int num, double ptL[][2],double pt[2]);
+void pull3(int num, double ptL[][3],double pt[3]);
+void pull4(int num, double ptL[][4],double pt[4]);
+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 prependd2(int from,int upto, double gdata[][2],double mat[][2]);
+int appendd1(double numvec[], double num);
+int appendi1(int numvec[],int num);
+int prependi1(int num, int numvec[]);
+int removei2(int nrmv,int mat[][2]);
+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]);
+void dispvec(int n,double dt[]);
+void dispveci(int n,int 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]);
+void dispmatd5(int from, int upto, double mat[][5]);
+void dispmatd6(int from, int upto, double mat[][6]);
+void dispmatd2all(double mat[][2]);
+void dispmatd3all(double mat[][3]);
+void dispmatd4all(double mat[][4]);
+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 readdata3(const char *fname, double out[][3]);
+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[]);
+int memberi(int a, int list[]);
+double dotprod(int dim, double p[], double q[]);
+double norm(int dim, double p[]);
+double dist(int dim, double p[], double q[]);
+void crossprod(int dim,double a[3],double b[3], double out[3]);
+void reflectpoint(double p1[2],double regard[][2],double p2[2]);
+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]);
+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]);
+void nearestptpt(double pa[2], double pL[][2], double ans[4]);
+int intersectselfPp(double data1[][2],double eps,double eps1s,double kL[][4]);
+int intersectcrvsPp(double data1[][2], double data2[][2],
+ double eps,double eps1s, double kL[][4]);
+int dropnumlistcrv(double qd[][2], double eps1, int out[]);
+int increasept2(double ptL[][2], int nn, double out[][2]);
+int increasept3(double ptL[][3], int nn, double out[][3]);
+void bezierpt(double t, double p0[2], double p3[2],
+ double p1[2], double p2[2], double ptout[2]);
+int bezier(double ptL[][2], double ctrL[][4], int num, double out[][2]);
+int ospline(int level, double ptL[][2], int num, double out[][2]);
+void bezierpt3(double t, double p0[3], double p3[3],
+ double p1[3], double p2[3], double ptout[3]);
+int bezier3(double ptL[][3], double ctrL[][6], int num, double out[][3]);
+int ospline3(int level, double ptL[][3], int num, double out[][3]);
+void segplot(double x1,double y1,double x2,double y2,double seg[][2]);
+void addsegplot(double x1,double y1,double seg[][2]);
+void circledata(double cx,double cy,double r, int num,double out[][2]);
+void intersectline(double p1[2],double v1[2], double p2[2],double v2[2],double out[4]);
+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]);
+
+const double Inf=99999;
+const int Infint=30000;
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h
new file mode 100644
index 00000000000..d9c0fbdc18c
--- /dev/null
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflib.h
@@ -0,0 +1,2098 @@
+// 17.06.16 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;
+ double tmpd[3];
+ FILE *fp;
+ fp=fopen(fname,"w");
+ 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++){
+ fprintf(fp,"start//\n");
+ fprintf(fp,"[");
+ ctr=0;
+ for(j=din[i][0]; j<=din[i][1]; j++){
+ pull3(j,out,tmpd);
+ fprintf(fp,"[ %7.5f, %7.5f, %7.5f]",tmpd[0],tmpd[1],tmpd[2]);
+ ctr++;
+ if(ctr==5){
+ fprintf(fp,"]//\n[");
+ ctr=0;
+ }else{
+ if(j<din[i][1]){
+ fprintf(fp,",");
+ }
+ }
+ }
+ fprintf(fp,"]//\n");
+ if(i<ndin){
+ fprintf(fp,"end//\n");
+ }
+ else{
+ fprintf(fp,"end////\n");
+ }
+ }
+ fclose(fp);
+ return 0;
+}
+
+int output3h(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");
+ if (fp == NULL) {
+ printf("cannot open\n");
+ return -1;
+ }
+ dataindexd3(3,out,din1);
+ out1[0][0]=0;
+ appendd3(0, din1[1][0],din1[1][1],out,out1);
+ dataindexd3(2, out1,din2);
+ fprintf(fp,"%s//\n",var);
+ for(i=1; i<=length2i(din2); i++){
+ fprintf(fp,"start//\n");
+ fprintf(fp,"[");
+ ctr=0;
+ for(j=din2[i][0]; j<=din2[i][1]; j++){
+ pull3(j,out1,tmpd);
+ fprintf(fp,"[ %7.5f, %7.5f, %7.5f]",tmpd[0],tmpd[1],tmpd[2]);
+ ctr++;
+ if(ctr==5){
+ if(j<din2[i][1]){
+ fprintf(fp,"]//\n[");
+ ctr=0;
+ }
+ }else{
+ if(j<din2[i][1]){
+ fprintf(fp,",");
+ }
+ }
+ }
+ fprintf(fp,"]//\n");
+ fprintf(fp,"end//\n");
+ }
+ out1[0][0]=0;
+ fprintf(fp,"%s//\n",varh); //17.06.16from
+ appendd3(0, din1[2][0],din1[2][1],out,out1);
+ dataindexd3(2, out1,din2);
+ for(i=1; i<=length2i(din2); i++){
+ fprintf(fp,"start//\n");
+ fprintf(fp,"[");
+ for(j=din2[i][0]; j<=din2[i][1]; j++){
+ pull3(j,out1,tmpd);
+ fprintf(fp,"[ %7.5f, %7.5f, %7.5f]",tmpd[0],tmpd[1],tmpd[2]);
+ ctr++;
+ if(ctr==5){
+ if(j<din2[i][1]){
+ fprintf(fp,"]//\n[");
+ ctr=0;
+ }
+ }else{
+ if(j<din2[i][1]){
+ fprintf(fp,",");
+ }
+ }
+ }
+ fprintf(fp,"]//\n");
+// if(i<ndin2){
+ fprintf(fp,"end//\n");
+// }
+// else{
+// fprintf(fp,"end////\n");
+// }
+ }
+ fprintf(fp,"//\n");
+ fclose(fp);
+ return 0;
+}
+
+int writedataC(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\n",out[i][0],out[i][1],out[i][2]);
+ }
+ 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;
+ FILE *fp;
+ fp=fopen(fname,"r");
+ if(fp==NULL){
+ printf("file not found\n");
+ return 1;
+ }
+ 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;
+ }
+ }
+ fclose(fp);
+ data[0][0]=ndata;
+ return ndata;
+}
+
+void parapt(double pt3[3],double pt2[2]){
+ if(fabs(pt3[0])>Inf-1){
+ pt2[0]=pt3[0]; pt2[1]=pt3[1];
+ }
+ else{
+ pt2[0]=-pt3[0]*sin(PHI)+pt3[1]*cos(PHI);
+ pt2[1]=-pt3[0]*cos(PHI)*cos(THETA)-pt3[1]*sin(PHI)*cos(THETA)+pt3[2]*sin(THETA);
+ }
+}
+
+int projpara(double p3[][3], double p2[][2]){
+ double pt3[3], pt2[2];
+ int i, nall, nall2;
+ nall=length3(p3);
+ p2[0][0]=0;
+ for(i=1; i<=nall; i++){
+ pull3(i,p3,pt3);
+ parapt(pt3,pt2);
+ nall2=addptd2(p2,pt2);
+ }
+ return nall2;
+}
+
+double invparapt(double ph, double fh[][2], double fk[][3],double out[3]){
+ int n, nfk=length3(fk), nph=length2(fh);
+ double s0=ph-n, ak[3], bk[3], ah[2], bh[2], pk[3];
+ if(nfk>2){
+ n=trunc(ph);
+ s0=ph-n;
+ if(ph>nph-Eps){
+ out[0]=fk[nfk][0]; out[1]=fk[nfk][1]; out[2]=fk[nfk][2];
+ return nph;
+ }
+ }
+ else{
+ n=1;
+ s0=ph-1;
+ }
+ ak[0]=fk[n][0]; ak[1]=fk[n][1]; ak[2]=fk[n][2];
+ bk[0]=fk[n+1][0]; bk[1]=fk[n+1][1]; bk[2]=fk[n+1][2];
+ ah[0]=fh[n][0]; ah[1]=fh[n][1];
+ bh[0]=fh[n+1][0]; bh[1]=fh[n+1][1];
+ out[0]=(1-s0)*ak[0]+s0*bk[0];
+ out[1]=(1-s0)*ak[1]+s0*bk[1];
+ out[2]=(1-s0)*ak[2]+s0*bk[2];
+ 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);
+ 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);
+ }
+ 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;
+ double eps=0.00001;
+ u1=u-eps/2; u2=u+eps/2;
+ if(u1<Urng[0]){u1=Urng[0];}
+ if(u2>Urng[1]){u2=Urng[1];}
+ surffun(ch,u1,v,pt1); parapt(pt1,pt2d1);
+ surffun(ch,u2,v,pt2); parapt(pt2,pt2d2);
+ dxu=(pt2d2[0]-pt2d1[0])/(u2-u1);
+ dyu=(pt2d2[1]-pt2d1[1])/(u2-u1);
+ v1=v-eps/2; v2=v+eps/2;
+ if(v1<Vrng[0]){v1=Vrng[0];}
+ if(v2>Vrng[1]){v2=Vrng[1];}
+ surffun(ch,u,v1,pt1); parapt(pt1,pt2d1);
+ surffun(ch,u,v2,pt2); parapt(pt2,pt2d2);
+ dxv=(pt2d2[0]-pt2d1[0])/(v2-v1);
+ dyv=(pt2d2[1]-pt2d1[1])/(v2-v1);
+ return dxu*dyv-dxv*dyu;
+};
+
+int envelopedata(short ch, double out[][2]){
+ int i, j, k, ctrq=0,ctr, nall;
+ 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;
+ ctr=0;
+ for (j = 0; j < Ndv; j++) {
+ vval1=Vrng[0]+j*dv;
+ vval2=Vrng[0]+(j+1)*dv;
+ uval1=Urng[0];
+ eval11=evlpfun(ch,uval1,vval1);
+ eval12=evlpfun(ch,uval1,vval2);
+ for (i = 0; i < Mdv; i++) {
+ uval2=Urng[0]+(i+1)*du;
+ eval21=evlpfun(ch,uval2,vval1);
+ eval22=evlpfun(ch,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;
+ 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++;
+ }
+ }
+ }
+ uval1=uval2;
+ eval11=eval21;
+ eval12=eval22;
+ 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 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;
+ int ntmp,ntmp1,ntmp2, ntmpnum, noutkL=0, noutk=0, nout=0;
+ int din[DsizeM][2], dind[DsizeM][2], drop[DsizeL], tmpnum[DsizeM];
+ double eps=pow(10.0,-6.0);
+ double vtmp, vtmp1,kaku,cva, ps[2], pe[2],v1[2],v2[2];
+ double tmp[2],tmp1[2],tmp2[2],pt3d[3],pt2d[2],q3d[3],q2d[2];
+ double ptk[DsizeM][3],pth[DsizeM][2], outk[DsizeM][3];
+ double ptxy[DsizeM][2],tmp3d[DsizeM][3];
+ double tmp3d1[DsizeM][3],tmp3d2[DsizeM][3],tmp2d[DsizeM][2];
+ outkL[0][0]=0;
+ cusp[0]=0;
+ Cuspsplitpt[0][0]=0;
+ nd=dataindexd2(0,pdata,din);
+ for(ng=1; ng<=nd; ng++){
+ ptk[0][0]=0; pth[0][0]=0;
+ st=din[ng][0]; en=din[ng][1];
+ ptxy[0][0]=0;
+ n=appendd2(0,st,en,pdata,ptxy);
+ npk=0; nph=0;
+ for(i=1; i<=en-st+1; i++){
+ pull2(st+i-1,pdata,tmp);
+ surffun(ch,tmp[0],tmp[1],pt3d);
+ parapt(pt3d,pt2d);
+ if(i==1){
+ npk=addptd3(ptk,pt3d);
+ nph=addptd2(pth,pt2d);
+ }
+ else{
+ pull2(nph,pth,tmp1);
+ if(dist(2,tmp1,pt2d)>eps){
+ npk=addptd3(ptk,pt3d);
+ nph=addptd2(pth,pt2d);
+ }
+ }
+ }
+ if(nph==0){
+ npk=0;
+ return 0;
+ };
+ pull2(1,pth,ps);
+ pull2(nph,pth,pe);
+ cflg=0;
+ if(dist(2,ps,pe)<Eps1){cflg=1;}
+ ncusp=0; cusp[0]=0;
+ cva=cos(10*M_PI/180);
+ for(i=1; i<=nph-1; i++){
+ if(nph==2){break;}
+ pull2(i,pth,tmp);
+ tmp[0]=pth[i][0]; tmp[1]=pth[i][1];
+ if(i==1){
+ if(cflg==0){continue;}
+ pull2(nph-1,pth,tmp1);
+ }
+ else{
+ pull2(i-1,pth,tmp1);
+ }
+ pull2(i+1,pth,tmp2);
+ v1[0]=tmp[0]-tmp1[0]; v1[1]=tmp[1]-tmp1[1];
+ v2[0]=tmp2[0]-tmp[0]; v2[1]=tmp2[1]-tmp[1];
+ vtmp=v1[0]*v2[0]+v1[1]*v2[1];
+ tmp[0]=0; tmp[1]=0;
+ vtmp=vtmp/(dist(2,v1,tmp)*dist(2,v2,tmp));
+ cuspflg=0;
+ if(vtmp<cva){
+ pt2d[0]=pth[i][0]; pt2d[1]=pth[i][1];
+ kaku=acos(vtmp)*180/M_PI;
+ if(v1[0]*v2[1]-v1[1]*v2[0]<0){
+ kaku=-kaku;
+ }
+ cuspflg=0;
+ for(j=i+1; j<=nph-1; j++){
+ if(fabs(kaku)>90){
+ cuspflg=1;
+ break;
+ }
+ pull2(j,pth,q2d);
+ if(dist(2,pt2d,q2d)>Eps1){
+ break;
+ }
+ v1[0]=q2d[0]-pth[j-1][0]; v1[1]=q2d[1]-pth[j-1][1];
+ v2[0]=pth[j+1][0]-q2d[0]; v2[1]=pth[j+1][1]-q2d[1];
+ vtmp=v1[0]*v2[0]+v1[1]*v2[1];
+ tmp[0]=0; tmp[1]=0;
+ vtmp=vtmp/(dist(2,v1,tmp)*dist(2,v2,tmp));
+ vtmp=acos(vtmp)*180/M_PI;
+ if(v1[0]*v2[1]-v1[1]*v2[0]<0){
+ vtmp=-vtmp;
+ }
+ kaku=kaku+vtmp;
+ }
+ if(cuspflg==1){
+ ntmp=trunc((i+j)*0.5);
+ i=j;
+ if(cusp[0]==0){
+ ncusp=appendi1(cusp,ntmp);
+ }
+ else{
+ ncusp=appendi1(cusp,ntmp);
+ }
+ }
+ }
+ }
+ if(cflg==0){
+ ncusp=prependi1(1,cusp);
+ ncusp=appendi1(cusp,nph);
+ }
+ else if(ncusp==0){
+ ncusp=prependi1(1,cusp);
+ ncusp=appendi1(cusp,nph);
+ }
+ else if(cusp[1]==1){
+ ncusp=appendi1(cusp,nph);
+ }
+ else{
+ ntmp=cusp[1];
+ ntmp1=nph;
+ tmp2d[0][0]=0;
+ n=appendd2(0,ntmp,ntmp1,pth,tmp2d);
+ n=appendd2(0,2,ntmp,pth,tmp2d);
+ pth[0][0]=0;
+ nph=appendd2(0,1,n,tmp2d,pth);
+ tmp3d[0][0]=0;
+ n=appendd3(0,ntmp,ntmp1,ptk,tmp3d);
+ n=appendd3(0,2,ntmp,ptk,tmp3d);
+ ptk[0][0]=0;
+ npk=appendd3(0,1,n,tmp3d,ptk);
+ for(i=1; i<=cusp[0]; i++){
+ cusp[i]=cusp[i]-ntmp+1;
+ }
+ ncusp=appendi1(cusp,nph);
+ }
+ if(ncusp==2){
+ pull2(nph,pth,tmp);
+ if(npk>=2){
+ ncusppt=addptd2(Cuspsplitpt,tmp);
+ noutkL=appendd3(2,1,npk,ptk,outkL);
+ }
+ continue;
+ }
+ outk[0][0]=0;
+ is=1;
+ for(i=1; i<=ncusp-1; i++){
+ ntmp1=cusp[is]; ntmp2=cusp[i+1];
+ pull2(ntmp1,pth,tmp1);
+ pull2(ntmp2,pth,tmp2);
+ if(dist(2,tmp1,tmp2)>Eps1){
+ tmp3d[0][0]=0;
+ n=appendd3(0,ntmp1,ntmp2,ptk,tmp3d);
+ noutk=appendd3(2,1,n,tmp3d,outk);
+ ncusppt=addptd2(Cuspsplitpt,tmp2);
+ is=i+1;
+ }
+ }
+ noutkL=appendd3(2,1,noutk,outk,outkL);
+ }
+ n=projpara(outkL,tmp2d);
+ ndrop=dropnumlistcrv(tmp2d,Eps1*0.5,drop);
+ nd=dataindexd3(2,outkL,din);
+ ndd=dataindexi1(drop,dind);
+ tmp3d2[0][0]=0; ntmp2=0;
+ for(i=1; i<=nd; i++){
+ tmp3d[0][0]=0; ntmp=0;
+ 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]);
+ }
+ tmp3d1[0][0]=0; ntmp1=0;
+ for(j=1; j<=ntmpnum; j++){
+ pt3d[0]=tmp3d[tmpnum[j]][0];
+ pt3d[1]=tmp3d[tmpnum[j]][1];
+ pt3d[2]=tmp3d[tmpnum[j]][2];
+ ntmp1=addptd3(tmp3d1,pt3d);
+ }
+ if(ntmp1>0){
+ ntmp2=appendd3(2,1,ntmp1,tmp3d1,tmp3d2);
+ }
+ }
+ outkL[0][0]=0; noutkL=0;
+ noutkL=appendd3(0,1,ntmp2,tmp3d2,outkL);
+ return noutkL;
+}
+
+int makexybdy(short ch,double ehL[][2]){
+ int i, nehL, ntmpd2;
+ double ptd3[3],ptd2[2],tmpd2[DsizeM][2],du,dv;
+ ehL[0][0]=0; nehL=0;
+ du=(Urng[1]-Urng[0])/Mdv;
+ dv=(Vrng[1]-Vrng[0])/Ndv;
+ if(DrawE==1){
+ tmpd2[0][0]=0; ntmpd2=0;
+ for(i=0; i<=Ndv; i++){
+ surffun(ch,Urng[1],Vrng[0]+i*dv,ptd3);
+ ntmpd2=push2(ptd3[0],ptd3[1],tmpd2);
+ }
+ nehL=appendd2(2,1,ntmpd2,tmpd2,ehL);
+ }
+ if(DrawN==1){
+ tmpd2[0][0]=0; ntmpd2=0;
+ for(i=0; i<=Mdv; i++){
+ surffun(ch,Urng[0]+i*du,Vrng[1],ptd3);
+ ntmpd2=push2(ptd3[0],ptd3[1],tmpd2);
+ }
+ nehL=appendd2(2,1,ntmpd2,tmpd2,ehL);
+ }
+ if(DrawW==1){
+ tmpd2[0][0]=0; ntmpd2=0;
+ for(i=0; i<=Ndv; i++){
+ surffun(ch,Urng[0],Vrng[0]+i*dv,ptd3);
+ ntmpd2=push2(ptd3[0],ptd3[1],tmpd2);
+ }
+ nehL=appendd2(2,1,ntmpd2,tmpd2,ehL);
+ }
+ if(DrawS==1){
+ tmpd2[0][0]=0; ntmpd2=0;
+ for(i=0; i<=Mdv; i++){
+ surffun(ch,Urng[0]+i*du,Vrng[0],ptd3);
+ ntmpd2=push2(ptd3[0],ptd3[1],tmpd2);
+ }
+ nehL=appendd2(2,1,ntmpd2,tmpd2,ehL);
+ }
+ return nehL;
+}
+
+void partitionseg(double fig[][2],double gL[][2], int ns, double parL[]){
+ int ii,jj,kk,n,din[DsizeM][2],din2[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){
+ dataindexd1(Otherpartition,din);
+ other[0]=0;
+ for(ii=din[ns][0]; ii<=din[ns][1]; ii++){
+ appendd1(other,Otherpartition[ii]);
+ }
+ }else{
+ other[0]=0;
+ }
+ parL[0]=0;
+ appendd1(parL,1);
+ appendd1(parL,length2(fig));
+ if(length1(other)>0){
+ for(ii=1; ii<=length1(other); ii++){
+ appendd1(parL,other[ii]);
+ }
+ }
+ dataindexd2(2,gL,din);
+ for(n=ns; n<=din[0][0]; n++){
+ g[0][0]=0;
+ appendd2(2,din[n][0],din[n][1],gL,g);
+ intersectcurvesPp(fig,g,20,kouL6);
+ kouL[0][0]=0;
+ for(ii=1;ii<=length6(kouL6);ii++){
+ pull6(ii,kouL6,tmpd6);
+ push4(tmpd6[0],tmpd6[1],tmpd6[2],tmpd6[3],kouL);
+ }
+ tmp1d1[0]=0;tmp2d1[0]=0;
+ for(jj=1;jj<=length4(kouL);jj++){
+ pull4(jj,kouL,tmpd4);
+ tmp=tmpd4[2];
+ if((tmp>1+Eps)&&(tmp<length2(fig)-Eps)){
+ appendd1(tmp1d1,tmp);
+ }
+ tmp=tmpd4[3];
+ if((tmp>1+Eps)&&(tmp<length2(g)-Eps)){
+ appendd1(tmp2d1,tmp);
+ }
+ }
+ for(jj=1;jj<=length1(tmp1d1);jj++){
+ appendd1(parL,tmp1d1[jj]);
+ }
+ 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]);
+ }
+ }
+ }
+ simplesort(parL);
+ tmp1d1[0]=0;
+ for(ii=1;ii<=length1(parL);ii++){
+ appendd1(tmp1d1,parL[ii]);
+ }
+ parL[0]=0;
+ tmp2=-1;
+ for(jj=1;jj<=length1(tmp1d1);jj++){
+ tmp=tmp1d1[jj];
+ if(fabs(tmp-tmp2)>Eps1){
+ appendd1(parL,tmp);
+ tmp2=tmp;
+ }
+ }
+}
+
+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);
+ if((fabs(Vec[1])>eps0)||(fabs(Vec[0])>eps0)){
+ out=(Vec[1]*(pt[0]-pa[0])-Vec[0]*(pt[1]-pa[1]));
+ }
+ else{
+ out=pt[0]-pa[0];
+ }
+ 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 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];
+// 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++){
+ 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;
+ 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)){
+ continue;
+ }
+ if((m1<=-eps0)&&(m2<=-eps0)){
+ continue;
+ }
+ if(fabs(m1)<eps0){
+ tmpd2[0]=p1[0]; tmpd2[1]=p1[1];
+ }
+ else if(fabs(m2)<eps0){
+ continue;
+ }
+ else{
+ tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
+ tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
+ }
+ nqL++;
+ qL[nqL][0]=tmpd2[0]; qL[nqL][1]=tmpd2[1];
+ }
+ if(nqL==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){
+ 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){
+ m1=pta[2]+v3/v2*(p1d3[1]-pta[1])-p1d3[2];
+ m2=pta[2]+v3/v2*(p2d3[1]-pta[1])-p2d3[2];
+ }
+ else{
+ m1=pta[1]-p1d3[1];
+ m2=pta[1]-p2d3[1];
+ }
+ if(m1*m2>=0){
+ if(((m1>eps0)&&(m2>eps0))||((m1<-eps0)&&(m2<-eps0))){
+ continue;
+ }
+ if(fabs(m1)<=eps0){
+ 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{
+ 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]);
+ }
+ tmpd3[0]=ptd3[0]-pta[0];
+ 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;
+ }
+ }
+ }
+ }
+ }
+ 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];
+ 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];
+ nfk=length3(fk);
+ nfh=projpara(fk, fh);
+ pull2(1,fh,p1);
+ pull2(nfh,fh,p2);
+ cspflg=1;
+ for(i=1; i<=ncusp; i++){
+ pull2(i, Cusppt, tmpd2);
+ if(dist(2,tmpd2,p1)<eps0){
+ if(cspflg==1){
+ cspflg=2;
+ }
+ else if(cspflg==3){
+ cspflg=6;
+ }
+ continue;
+ }
+ if(dist(2,tmpd2,p2)<eps0){
+ if(cspflg==1){
+ cspflg=3;
+ }
+ else if(cspflg==2){
+ cspflg=6;
+ }
+ 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);
+ 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);
+ }
+ }
+ figL[0][0]=0; nfigL=0;
+ figkL[0][0]=0; nfigkL=0;
+ for(i=1; i<=nseL; i++){
+ pointoncurve(paL[seL[i]], fh, Eps, tmp1d2);
+ pointoncurve(paL[seL[i]+1], fh, Eps, tmp2d2);
+ if(i==1){
+ copyd(0,1,tmp1d2,pd2);
+ sp=paL[seL[i]];
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ else{
+ if(memberi(seL[i]-1,seL)==1){
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ else{
+ ntmp=partcrv(sp,sq,fh,tmpmd2);
+ nfigL=appendd2(2,1,ntmp,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);
+ copyd(0,1,tmp1d2,pd2);
+ sp=paL[seL[i]];
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ }
+ }
+ if(nseL>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);
+ 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);
+ }
+ else{
+ nfigL=appendd2(2,1,nfh, fh, figL);
+ nfigkL=appendd3(2,1,nfk,fk,figkL);
+ }
+ }
+ tmpi1[0]=0; ntmp=0;
+ for(i=1; i<=npaL-1; i++){
+ if(memberi(i,seL)==0){
+ ntmp=appendi1(tmpi1,i);
+ }
+ }
+ copyi(0,ntmp,tmpi1,seL);
+ seL[0]=ntmp; nseL=ntmp;
+ Hiddendata[0][0]=0; nhidden=0;
+ for(i=1; i<=nseL; i++){
+ dtmp=paL[seL[i]];
+ pointoncurve(dtmp,fh,Eps,tmp1d2);
+ dtmp=paL[seL[i]+1];
+ pointoncurve(dtmp,fh,Eps, tmp2d2);
+ if(i==1){
+ copyd(0,1,tmp1d2,pd2);
+ sp=paL[seL[i]];
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ else{
+ if(memberi(seL[i]-1,seL)==1){
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ 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);
+ copyd(0,1,tmp1d2,pd2);
+ sp=paL[seL[i]];
+ copyd(0,1,tmp2d2,qd2);
+ sq=paL[seL[i]+1];
+ }
+ }
+ }
+ if(nseL>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);
+ }
+ else{
+ nhidden=appendd3(2,1,nfk,fk,Hiddendata);
+ }
+ }
+ return nfigkL;
+}
+
+int 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;
+ 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;
+ if(DrawE==1){
+ tmpd2[0][0]=0; nall=0;
+ for(j=0; j<=Ndv; j++){
+ nall=add2(tmpd2, Urng[1], Vrng[0]+j*dv);
+ }
+ nall=surfcurve(ch,2, tmpd2,tmpd3);
+ nekL=appendd3(2,1,nall,tmpd3,ekL);
+ }
+ if(DrawN==1){
+ tmpd2[0][0]=0; nall=0;
+ for(j=0; j<=Mdv; j++){
+ nall=add2(tmpd2, Urng[0]+j*du, Vrng[1]);
+ }
+ nall=surfcurve(ch,2, tmpd2,tmpd3);
+ nekL=appendd3(2,1,nall,tmpd3,ekL);
+ }
+ if(DrawW==1){
+ tmpd2[0][0]=0; nall=0;
+ for(j=0; j<=Ndv; j++){
+ nall=add2(tmpd2, Urng[0], Vrng[0]+j*dv);
+ }
+ nall=surfcurve(ch,2, tmpd2,tmpd3);
+ nekL=appendd3(2,1,nall,tmpd3,ekL);
+ }
+ if(DrawS==1){
+ tmpd2[0][0]=0; nall=0;
+ for(j=0; j<=Mdv; j++){
+ nall=add2(tmpd2, Urng[0]+j*du, Vrng[0]);
+ }
+ nall=surfcurve(ch,2, tmpd2,tmpd3);
+ nekL=appendd3(2,1,nall,tmpd3,ekL);
+ }
+ if(nekL>0){
+ nfkL=appendd3(2,1,nekL,ekL,fkL);
+ }
+ nfall=projpara(fkL,fall);
+ nfbdxy=makexybdy(ch,fbdxy);
+ ndin=dataindexd2(2,fall,din);
+ Borderpt[0][0]=0;
+ Otherpartition[0]=0;
+ for(i=1; i<=ndin-1; i++){
+ Otherpartition[i]=Inf; Otherpartition[0]++;
+ }
+ Borderhiddendata[0][0]=0;
+ ndin2=dataindexd3(2,fkL,din2);
+ for(i=1; i<= ndin2; 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);
+ if(length1(par)>2){
+ ntmp=npar-1;
+ tmpd2[0][0]=0; ntmpd2=0;
+ for(j=2; j<=ntmp; j++){
+ pull2(j,Partitionpt,p);
+ ntmpd2=addptd2(tmpd2, p);
+ }
+ ntmp=appendd2(0, 1,ntmpd2,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);
+ }
+ ntmp=length3(Hiddendata);
+ if(ntmp>0){
+ appendd3(2,1,ntmp,Hiddendata,Borderhiddendata);
+ }
+ }
+ return nfsL;
+}
+
+int dropnumlistcrv3(double qd[][3], double eps1, int out[]){
+ int i,j,k,nout,nall=length3(qd), nd, se, en, npd, nptL;
+ int din[DsizeM][2],ptL[DsizeL];
+ double eps=pow(10.0,-6.0), pd[DsizeL][3], p[3], tmp2d[3];
+ out[0]=0; nout=0;
+ nd=dataindexd3(2,qd,din);
+ for(j=1; j<=nd; j++){
+ se=din[j][0]; en=din[j][1];
+ pd[0][0]=0; npd=0;
+ npd=appendd3(0,se,en,qd,pd);
+ ptL[0]=0; nptL=0;
+ nptL=appendi1(ptL,1);
+ 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);
+ 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);
+ }
+// if(nptL==1){ptL[0]=0; nptL=0;}//180318db
+ if(nout>0){
+ nout=appendi1(out,Infint);
+ }
+ for(i=1; i<=nptL; i++){
+ nout=appendi1(out,ptL[i]);
+ }
+ }
+ out[0]=nout;
+ return nout;
+}
+
+int 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;
+ pts[0][0]=0;
+ dataindexd2(2,tmp3md,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);
+ addptd2(tmp1md,tmpd);
+ surffun(ch,tmpd[0],tmpd[1],tmpd3);
+ addptd3(tmp2md3,tmpd3);
+ }
+ dropnumlistcrv3(tmp2md3,Eps1,nlist);
+ if(nlist[0]==1){
+ pull2(nlist[1],tmp1md,tmpd);
+ addptd2(pts,tmpd);
+ }else{
+ tmpmd[0][0]=0;
+ for(kk=1;kk<=nlist[0];kk++){
+ pull2(nlist[kk],tmp1md,tmpd);
+ addptd2(tmpmd,tmpd);
+ }
+ if(length2(tmpmd)>0){
+ appendd2(2,1,length2(tmpmd),tmpmd,out3md);
+ }
+ }
+ }
+ tmp3md[0][0]=0;
+ appendd2(2,1,length2(pts),pts,tmp3md);
+ pts[0][0]=0;
+ for(ii=1;ii<=length2(tmp3md);ii++){
+ pull2(ii,tmp3md,tmpd);
+ surffun(ch,tmpd[0],tmpd[1],tmp1d3);
+ flg=0;
+ for(jj=1;jj<=length2(pts);jj++){
+ pull2(jj,pts,tmpd);
+ surffun(ch,tmpd[0],tmpd[1],tmp2d3);
+ if(dist(3,tmp1d3,tmp2d3)<Eps1){
+ flg=1; break;
+ }
+ }
+ if(flg==0){
+ pull2(ii,tmp3md,tmpd);
+ addptd2(pts,tmpd);
+ }
+ }
+ 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);
+ 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);
+ eps=1;
+ for(kk=1;kk<=length2(tmp2md)-1;kk++){
+ pull2(kk,tmp2md,tmpd);
+ surffun(ch,tmpd[0],tmpd[1],tmp1d3);
+ pull2(kk+1,tmp2md,tmpd);
+ surffun(ch,tmpd[0],tmpd[1],tmpd3);
+ if(eps>dist(3,tmp1d3,tmpd3)){
+ eps=dist(3,tmp1d3,tmpd3);
+ }
+ }
+ 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);
+ }
+ 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);
+ 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);
+ 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]));
+ }
+ else{
+ Zval[j][i]=pt[0]-pa[0];
+ }
+ }
+ }
+ 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 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;
+ 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)){
+ continue;
+ }
+ if((m1<=-eps0)&&(m2<=-eps0)){
+ continue;
+ }
+ if(fabs(m1)<eps0){
+ copyd(0,1,p1d2,tmpd2);
+ }
+ else if(fabs(m2)<eps0){
+ continue;
+ }
+ else{
+ tmpd2[0]=1/(m1-m2)*(-m2*p1d2[0]+m1*p2d2[0]);
+ tmpd2[1]=1/(m1-m2)*(-m2*p1d2[1]+m1*p2d2[1]);
+ }
+ nqL=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);
+ v1=vec[0]; v2=vec[1]; v3=vec[2];
+ if(fabs(v1)>eps0){
+ 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){
+ m1=pta[2]+v3/v2*(p1d3[1]-pta[1])-p1d3[2];
+ m2=pta[2]+v3/v2*(p2d3[1]-pta[1])-p2d3[2];
+ }
+ else{
+ m1=pta[1]-p1d3[1];
+ m2=pta[1]-p2d3[1];
+ }
+ if(m1*m2>=-eps0){
+ if(((m1>eps0)&&(m2>eps0)) || ((m1<-eps0)&&(m2<-eps0))){
+ 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);
+ }
+ }
+ 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]);
+ }
+ 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);
+ }
+ }
+ }
+ }
+ 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);
+ }
+ }
+ 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];
+ double fh[DsizeM][2], parL[DsizeM], tmpmd3[DsizeM][3];
+ double ptL[DsizeM][3], tmpd2[2], tmpd3[3], tmp1d1, tmp2d1, dt;
+ 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);
+ 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++){
+ 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);
+ parapt(tmpd3,tmpd2);
+ tmp1d1=paramoncurve(tmpd2,i,fh);
+ tmp2d1=Inf;
+ for(k=1; k<=npar; 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);
+ }
+ }
+ simplesort(parL);
+ }
+ }
+ nohiddenpara2(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;
+}
+
+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);
+ }
+ }
+ else{
+ for(i=1; i<=n; i++){
+ nu=appendd1(udiv,udv[i]);
+ }
+ }
+ 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);
+ }
+ }
+ 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);
+ }
+ }
+ nohiddenpara2(par,fk,1, fs);
+ appendd3(2,1,length3(Hiddendata), Hiddendata,fsLh);
+ appendd3(2,1,length3(fs),fs,fsL);
+ }
+ }
+ push3(Inf,3,0,fsL);
+ nfsL=appendd3(0,1,length3(fsLh),fsLh,fsL);
+ push3(Inf,3,0,fsL);
+ return nfsL;
+}
+
+/*
+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 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);
+ }
+ }
+ 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 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;
+ int i, n, nall;
+ out[0][0]=0;
+ n=length3(curve);
+ for(i=1; i<=n; i++){
+ pull3(i, curve, pt);
+ if(pt[0]>Inf-1){
+ xz=pt[0]; yz=pt[1]; zz=pt[2];
+ }
+ else{
+ xz=-pt[0]*SP+pt[1]*CP;
+ yz=-pt[0]*CP*CT-pt[1]*SP*CT+pt[2]*ST;
+ zz=pt[0]*CP*ST+pt[1]*SP*ST+pt[2]*CT;
+ }
+ nall=push3(xz,yz,zz,out);
+ }
+ return nall;
+}
+
+int kukannozoku(double jokyo[2], double kukanL[][2], double res[][2]){
+ double t1, t2, tmpd2[2], ku[2];
+ int i, j, ii, nn, flg, nres;
+ nn=length2(kukanL);
+ t1=jokyo[0];
+ t2=jokyo[1];
+ pull2(1,kukanL,tmpd2);
+ t1=fmax(t1,tmpd2[0]);
+ pull2(nn,kukanL,tmpd2);
+ t2=fmin(t2,tmpd2[1]);
+ res[0][0]=0; nres=0;
+ flg=0;
+ for(ii=1; ii<=nn; ii++){
+ pull2(ii,kukanL,ku);
+ if(flg==0){
+ if(ku[1]<t1){
+ nres=addptd2(res,ku);
+ }
+ else{
+ flg=1;
+ if(ku[0]<t1-Eps){
+ tmpd2[0]=ku[0]; tmpd2[1]=t1;
+ nres=addptd2(res,tmpd2);
+ }
+ if(ku[1]>t2+Eps){
+ tmpd2[0]=t2; tmpd2[1]=ku[1];
+ nres=addptd2(res,tmpd2);
+ }
+ }
+ }
+ else if(flg==1){
+ if(ku[1]<t2){
+ continue;
+ }
+ else{
+ flg=2;
+ if(ku[0]<t2-Eps){
+ ku[0]=t2;
+ }
+ nres=addptd2(res,ku);
+ }
+ }
+ else{
+ nres=addptd2(res,ku);
+ }
+ }
+ return nres;
+}
+
+int skeletondata3(double data[][3], double r00,
+ double eps1, double eps2, double allres[][3]){
+ double objL[DsizeLL][3], pltL[DsizeLL][3];
+ double clipL[DsizeL][5],koc2d6[DsizeM][6],tmpd6[6];
+ double pts1[DsizeL][2],koc[DsizeM][4], koc2[DsizeM][4], kukanL[DsizeL][2];
+ double r0, t1, t2, z1, z2, tt, rr, origin[2], te, hh, ku[2], za, zb, ds[DsizeM],dmin;
+ double pt[2], pta[2], ptb[2], ptq[2];
+ double ptL[DsizeL][3], ptL2[DsizeL][2], res2[DsizeL][2], allres2[DsizeL][2];
+ double obj[DsizeL][3], plt[DsizeL][3];
+ double objp[DsizeL][3], pltp[DsizeL][3];
+ double obj2d[DsizeL][2], plt2d[DsizeL][2];
+ double tmpd1, tmp2d1, tmp3d1, tmp1d2[2], tmp2d2[2],tmpd3[3];
+ double tmpmd2[DsizeL][2], tmpmd3[DsizeL][3], tmpd4[4], tmp1d4[4];
+ int j, ii, jj, nobj, nplt, nall, din1[DsizeM][2], din2[DsizeM][2];
+ int nkoc, nkoc2, nn,nn1,nn2, nclipL, flg, ntmp;
+ char str1[100], str2[100];
+ objL[0][0]=0; pltL[0][0]=0;
+ dataindexd3(3,data,din1);
+ appendd3(0, din1[1][0],din1[1][1],data,objL);
+ appendd3(0, din1[2][0],din1[2][1],data,pltL);
+ allres[0][0]=0;
+ allres2[0][0]=0;
+ dataindexd3(2,objL,din1);
+ dataindexd3(2,pltL,din2);
+ for(nobj=1; nobj<=length2i(din1); nobj++){
+ tmpmd3[0][0]=0;
+ ntmp=appendd3(2,din1[nobj][0],din1[nobj][1],objL,tmpmd3);
+ nn=ceil((fmax(Mdv,Ndv)*2-1)/(ntmp-1));
+ increasept3(tmpmd3,nn,obj);
+ appendd3(0,ntmp,ntmp,tmpmd3,obj);
+ projcoordcurve(obj,objp);
+ obj2d[0][0]=0;
+ for(j=1; j<=length3(objp); j++){
+ push2(objp[j][0],objp[j][1],obj2d);
+ }
+ clipL[0][0]=0; nclipL=0;
+ for(nplt=1; nplt<=length2i(din2); nplt++){
+ tmpmd3[0][0]=0;
+ ntmp=appendd3(2,din2[nplt][0],din2[nplt][1],pltL,tmpmd3);
+ nn=ceil((fmax(Mdv,Ndv)*2-1)/(ntmp-1));
+ increasept3(tmpmd3,nn,plt);
+ appendd3(0,ntmp,ntmp,tmpmd3,plt);
+ projcoordcurve(plt,pltp);
+ plt2d[0][0]=0;
+ for(j=1; j<=length3(pltp); j++){
+ push2(pltp[j][0],pltp[j][1],plt2d);
+ }
+ intersectcurvesPp(obj2d,plt2d,20,koc2d6); //180324from
+ nkoc2=length6(koc2d6);
+ koc2[0][0]=0;
+ for(ii=1;ii<=nkoc2;ii++){
+ pull6(ii,koc2d6,tmpd6);
+ push4(tmpd6[0],tmpd6[1],tmpd6[2],tmpd6[3],koc2);
+ }//180324to
+ if(nkoc2>0){
+ koc[0][0]=0; nkoc=0;
+ pull2(1,plt2d,pta);
+ pull2(length2(plt2d),plt2d,ptb);
+ for(j=1; j<=nkoc2; j++){
+ pt[0]=koc2[j][0]; pt[1]=koc2[j][1];
+ nearestptpt(pt, plt2d, tmpd4);
+ koc2[j][3]=tmpd4[2];
+ pt[0]=tmpd4[0]; pt[1]=tmpd4[1];
+ nearestptpt(pt, obj2d, tmpd4);
+ koc2[j][0]=tmpd4[0]; koc2[j][1]=tmpd4[1];
+ koc2[j][2]=tmpd4[2];
+ ds[j]=tmpd4[3];
+ }
+ pull4(1,koc2,tmp1d4);
+ dmin=ds[1];
+ for(j=2; j<=nkoc2; j++){
+ if(fabs(koc2[j][2]-tmp1d4[2])<4){
+ if(ds[j]<dmin){
+ pull4(j,koc2,tmp1d4);
+ dmin=ds[j];
+ }
+ }
+ else{
+ addptd4(koc,tmp1d4);
+ pull4(j,koc2,tmp1d4);
+ dmin=ds[j];
+ }
+ }
+ nkoc=addptd4(koc,tmp1d4);
+/*
+ if(tmpd4[3]<2*eps1){
+ if((dist(2,pta,pt)<eps2)&&(tmpd4[2]<0.1*length2(plt2d))){
+ continue;
+ }
+ if((dist(2,ptb,pt)<eps2)&&(tmpd4[2]>0.9*length2(plt2d))){
+ continue;
+ }
+ pull4(j,koc2,tmpd4);
+ nkoc=addptd4(koc,tmpd4);
+ }
+ }
+*/
+ for(jj=1; jj<=nkoc; jj++){
+ pt[0]=koc[jj][0]; pt[1]=koc[jj][1];
+ pta[0]=plt2d[1][0]; pta[1]=plt2d[1][1];
+ ptb[0]=plt2d[length2(plt2d)][0]; ptb[1]=plt2d[length2(plt2d)][1];
+ if((dist(2,pt,pta)<eps2)&&(koc[jj][3]<0.3*length2(plt2d))){
+ continue;
+ }
+ if((dist(2,pt,ptb)<eps2)&&(koc[jj][3]>0.7*length2(plt2d))){
+ continue;
+ }
+ nn1=floor(koc[jj][2]);
+ nn2=koc[jj][3];
+ pta[0]=objp[nn1][0]; pta[1]=objp[nn1][1];
+ za=objp[nn1][2];
+ ptb[0]=objp[nn1+1][0]; ptb[1]=objp[nn1+1][1];
+ zb=objp[nn1+1][2];
+ if(dist(2,pta,ptb)<eps1){
+ continue;
+ }
+ t1=dist(2,pta,pt)/dist(2,pta,ptb);
+ z1=(1-t1)*za+t1*zb;
+ pta[0]=pltp[nn2][0]; pta[1]=pltp[nn2][1];
+ za=tmpmd3[nn2][2];
+ ptb[0]=pltp[nn2+1][0]; ptb[1]=pltp[nn2+1][1];
+ zb=tmpmd3[nn2+1][2];
+ if(dist(2,pta,ptb)<Eps){
+ continue;
+ }
+ t2=dist(2,pta,pt)/dist(2,pta,ptb);
+ z2=(1-t2)*za+t2*zb;
+ if(z1<z2-eps1){
+ if(nclipL==0){
+ tmpd1=1;
+ }
+ }
+ else{
+ tmpd1=Inf;
+ for(j=1; j<=nclipL; j++){
+ tmp2d1=pow(clipL[j][0]-pt[0],2.0);
+ tmp2d1=tmp2d1+pow(clipL[j][1]-pt[1],2.0);
+ tmpd1=fmin(tmpd1,tmp2d1);
+ }
+ }
+ if(tmpd1>pow(Eps,2.0)){
+ pts1[0][0]=0;
+ pointoncurve(nn1,obj2d,Eps,tmp1d2);
+ addptd2(pts1,tmp1d2);
+ pointoncurve(nn1+1,obj2d,Eps,tmp1d2);
+ addptd2(pts1,tmp1d2);
+ tmp1d2[0]=pts1[2][0]-pts1[1][0];
+ tmp1d2[1]=pts1[2][1]-pts1[1][1];
+ tmp2d2[0]=ptb[0]-pta[0];
+ tmp2d2[1]=ptb[1]-pta[1];
+ tmp3d1=dotprod(2,tmp1d2,tmp2d2);
+ origin[0]=0; origin[1]=0;
+ tmp3d1=tmp3d1/dist(2,tmp1d2,origin)/dist(2,tmp2d2,origin);
+ tmpd1=1-0.5*pow(tmp3d1,2.0);
+ r0=0.075*r00;
+ nclipL=push5(pt[0],pt[1], nn1,t1,r0/tmpd1, clipL);
+ }
+ }
+ }
+ }
+ kukanL[0][0]=0;
+ te=length3(obj);
+ push2(1.0,te, kukanL);
+ for(ii=1; ii<=nclipL; ii++){
+ pt[0]=clipL[ii][0]; pt[1]=clipL[ii][1];
+ nn=clipL[ii][2];
+ tt=nn+clipL[ii][3];
+ rr=clipL[ii][4];
+ flg=0;
+ for(jj=nn; jj>=1; jj--){
+ pointoncurve(jj,obj2d,Eps,ptq);
+ if(dist(2,pt,ptq)<rr){
+ continue;
+ }
+ flg=jj;
+ break;
+ }
+ if(flg==0){
+ t1=1;
+ }
+ else{
+ t1=flg; t2=tt;
+ hh=t2-t1;
+ for(j=1; j<=10; j++){
+ hh=hh*0.5;
+ pointoncurve(t1+hh,obj2d,Eps,ptq);
+ if(dist(2,pt,ptq)<rr){
+ t2=t2-hh;
+ }
+ else{
+ t1=t1+hh;
+ }
+ }
+ }
+ ku[0]=t1;
+ flg=0;
+ for(jj=nn+1; jj<=te; jj++){
+ pointoncurve(jj,obj2d,Eps,ptq);
+ if(dist(2,pt,ptq)<rr){
+ continue;
+ }
+ flg=jj;
+ break;
+ }
+ if(flg==0){
+ t2=te;
+ }
+ else{
+ t1=tt; t2=flg;
+ hh=t2-t1;
+ for(j=1; j<=10; j++){
+ hh=hh*0.5;
+ pointoncurve(t1+hh,obj2d,Eps,ptq);
+ if(dist(2,pt,ptq)<rr){
+ t1=t1+hh;
+ }
+ else{
+ t2=t2-hh;
+ }
+ }
+ }
+ ku[1]=t2;
+ ntmp=kukannozoku(ku,kukanL,tmpmd2);
+ kukanL[0][0]=0;
+ appendd2(0,1,ntmp,tmpmd2,kukanL);
+ }
+ for(jj=1; jj<=length2(kukanL); jj++){
+ tmp1d2[0]=kukanL[jj][0]; tmp1d2[1]=kukanL[jj][1];
+ t1=kukanL[jj][0]; nn1=floor(t1);
+ t2=kukanL[jj][1]; nn2=floor(t2);
+ ptL2[0][0]=0;
+ ptL[0][0]=0;
+ if(t1-nn1<1-Eps){
+ pointoncurve(t1,obj2d,Eps,tmp1d2);
+ addptd2(ptL2,tmp1d2);
+ invparapt(t1, obj2d, obj, tmpd3);
+ addptd3(ptL,tmpd3);
+ }
+ for(j=nn1+1; j<=nn2; j++){
+ pointoncurve(j,obj2d,Eps,tmp1d2);
+ addptd2(ptL2,tmp1d2);
+ invparapt(j, obj2d, obj, tmpd3);
+ addptd3(ptL,tmpd3);
+ }
+ if(t2-nn2>Eps){
+ pointoncurve(t2,obj2d,Eps,tmp1d2);
+ addptd2(ptL2,tmp1d2);
+ invparapt(t2, obj2d, obj, tmpd3);
+ addptd3(ptL,tmpd3);
+ }
+ if(length2(ptL2)>1){
+ appendd2(2,1,length2(ptL2),ptL2,allres2);
+ appendd3(2,1,length3(ptL),ptL,allres);
+ }
+ }
+ }
+ 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);
+}
+*/
+
+
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h
new file mode 100644
index 00000000000..df31a4d812f
--- /dev/null
+++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/surflibhead.h
@@ -0,0 +1,58 @@
+double Xval[DsizeS], Yval[DsizeS], Zval[DsizeS][DsizeS];
+
+//////////////////////
+short Ch=1;
+/////////////////////
+
+//double Implicitdata[DsizeL][2]; //180320
+double Cuspsplitpt[DsizeL][2];
+double Cusppt[DsizeM][2];
+double Otherpartition [DsizeM];
+double Partitionpt[DsizeM][2];
+double Cuspdata[DsizeM];
+double Wirept[DsizeM][3];
+double Hiddendata[DsizeL][3];
+double Borderpt[DsizeM][2];
+double Borderhiddendata[DsizeL][3];
+double Surf[DsizeL][3];
+double Crvdata[DsizeL][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 writedataC(const char *fname, double out[][3]);
+int readdataC(const char *fname, double data[][3]);
+void parapt(double pt3[3],double pt2[2]);
+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]);
+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]);
+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]);
+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);