// 180430
//    crv2sfparadata,wireparadata added
//    output3h,output3 chanded
// 180426 readdataC,writedataC changed
// 180420 newly rewrited
// 170616 crvsfparadata changed ( num of devision )

int output3(short head, const char *wa,const char *var, const char *fname,double out[][3]){
 int din[DsizeS][2],i,j,ctr;
 double tmpd[3];
 FILE *fp;
 fp=fopen(fname,wa);
 if (fp == NULL) {
   printf("cannot open\n");
   return -1;
 }
 dataindexd3(2, out,din);
 if(head==1){
   fprintf(fp,"%s//\n",var);
 }
 for(i=1; i<=length2i(din); i++){
   fprintf(fp,"start//\n");
   fprintf(fp,"[");
   ctr=0;
   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;
       if(j<din[i][1]){
         fprintf(fp,"[");
       }
     }else{
       if(j<din[i][1]){
         fprintf(fp,",");
       }else{
         fprintf(fp,"]//\n");
       }
     }
   }
   fprintf(fp,"end//\n");
 }
 fclose(fp);
 return 0;
}

int output3h(const char *wa,const char *var, const char *varh, const char *fname, double out[][3]){
 // 17.06.09
 double out1[DsizeLL][3], tmpd[3];
 int din1[DsizeM][2],din2[DsizeM][2], i, j,ctr;
 FILE *fp;
 fp=fopen(fname,wa);
 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
 if(din1[0][0]>1){ //180613from
   appendd3(0, din1[2][0],din1[2][1],out,out1);
 }//180613to
 dataindexd3(2, out1,din2);
 for(i=1; i<=length2i(din2); i++){
   fprintf(fp,"start//\n");
   fprintf(fp,"[");
   ctr=0; //181105
   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");
 }
//  fprintf(fp,"//\n"); //180531
 fclose(fp);
 return 0;
}

void outputend(const char *fname){
 FILE *fp;
 fp=fopen(fname,"a");
 fprintf(fp,"//\n");
 fclose(fp);
}

int writedataC(const char *fname, double out[][3]){
 int i, nall;
 FILE *fp;
 nall=length3(out);
 fp=fopen(fname,"w");
 for(i=1; i<=nall; i++){
   fprintf(fp,"%7.5f %7.5f %7.5f %6s",out[i][0],out[i][1],out[i][2],"-99999");
   if(i<nall){
     fprintf(fp,"\n");
   };
 }
 fclose(fp);
 return nall;
}

int readdataC(const char *fname, double data[][3]){
 int i, j, ndata=0;
 double tmpd4[4];
 FILE *fp;
 fp=fopen(fname,"r");
 if(fp==NULL){
   printf("file not found\n");
   return 1;
 }
 data[0][0]=0;
 while( ! feof(fp) && ndata<DsizeL){
   for(j=0;j<=3;j++){
     fscanf(fp,"%lf",&tmpd4[j]);
   }
   push3(tmpd4[0],tmpd4[1],tmpd4[2],data);
 }
 fclose(fp);
 return length3(data);
}

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);
 }
}

double zparapt(double p[3]){
 double x=p[0], y=p[1], z=p[2], out;
 out=x*cos(PHI)*sin(THETA)+y*sin(PHI)*sin(THETA)+z*cos(THETA);
 return out;
}

int projpara(double p3[][3], double p2[][2]){
 double pt3[3], pt2[2];
 int i, nall, nall2;
 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;
}

void surfcurve(short ch,double crv[][2], double pdata[][3]){
 double pd2[2], pt[3];
 int i;
 pdata[0][0]=0;
 for(i=1; i<=length2(crv); i++){
   pull2(i,crv,pd2);
   surffun(ch,pd2[0],pd2[1],pt);
   addptd3(pdata,pt);
 }
}

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])/Ndv; //180427
 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 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(ntmp,cusp);
         }
         else{
           ncusp=appendi1(ntmp,cusp);
         }
       }
     }
   }
   if(cflg==0){
     ncusp=prependi1(1,cusp);
     ncusp=appendi1(nph,cusp);
   }
   else if(ncusp==0){
      ncusp=prependi1(1,cusp);
      ncusp=appendi1(nph,cusp);
   }
   else if(cusp[1]==1){
     ncusp=appendi1(nph,cusp);
   }
   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(nph,cusp);
   }
   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(drop[j],tmpnum);
   }
   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];
 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>0){
   dataindexd1(Otherpartition,din);
   other[0]=0;
   for(ii=din[ns][0]; ii<=din[ns][1]; ii++){
     appendd1(Otherpartition[ii],other);
   }
 }else{
   other[0]=0;
 }
 parL[0]=0;
 appendd1(1,parL);
 appendd1(length2(fig),parL);
 if(length1(other)>0){
   for(ii=1; ii<=length1(other); ii++){
     appendd1(other[ii],parL);
   }
 }
 dataindexd2(2,gL,din);
 for(n=fmax(ns,1); n<=din[0][0]; n++){
   g[0][0]=0;
   appendd2(2,din[n][0],din[n][1],gL,g);
       intersectcurvesPp(fig,g,20,kouL6);
   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(tmp,tmp1d1);
     }
     tmp=tmpd4[3];
     if((tmp>1+Eps)&&(tmp<length2(g)-Eps)){
       appendd1(tmp,tmp2d1);
     }
   }
   for(jj=1;jj<=length1(tmp1d1);jj++){
     appendd1(tmp1d1[jj],parL);
   }
   if((ns>0)&&(n>ns)&&(length1(tmp2d1)>0)){
     dataindexd1(Otherpartition,din);
     tmp1d1[0]=0;
     for(jj=din[n][0];jj<=din[n][1];jj++){
       appendd1(Otherpartition[jj],tmp1d1);
     }
     replacelistd1(n,Otherpartition,tmp2d1);
   }
 }
 simplesort(parL);
 tmp1d1[0]=0;
 for(ii=1;ii<=length1(parL);ii++){
   appendd1(parL[ii],tmp1d1);
 }
 parL[0]=0;
 tmp2=-1;
 for(jj=1;jj<=length1(tmp1d1);jj++){
   tmp=tmp1d1[jj];
   if(fabs(tmp-tmp2)>Eps1){
     appendd1(tmp,parL);
     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;
}

int pthiddenQ(short ch,double pta[3], int uveq,double out[4]){
 int i, j, k;
 double Vec[3]={sin(THETA)*cos(PHI),sin(THETA)*sin(PHI),cos(THETA)};
 double v1,v2,v3;
 double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
 double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
 double pl[5][2], vl[5], ql[11][2],p1[2],p2[2],m1,m2;
 double ptuv[2], puv1[2], puv2[2],p1d3[3], p2d3[3], ptd3[3];
 double tmpd1, tmpd2[2],tmpd3[3],tmp1d3[3];
//     Out=1 : hidden
//     SL : List of segments near PtA+Vec
 out[3]=-Inf;
 for(j=0; j<=Ndv-1; j++){
   vval1=Vrng[0]+j*dv;
   vval2=Vrng[0]+(j+1)*dv;
   uval1=Urng[0];
   eval11=funpthiddenQ(ch,uval1,vval1,pta);
   eval12=funpthiddenQ(ch,uval1,vval2,pta);
   for(i=0; i<=Mdv-1; i++){
     uval2=Urng[0]+(i+1)*du;
     eval21=funpthiddenQ(ch,uval2,vval1,pta);
     eval22=funpthiddenQ(ch,uval2,vval2,pta);
         pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
     pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
     pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
     pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
         pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
     ql[0][0]=0;
     for(k=0; k<=3; k++){
       pull2(k,pl,p1);
       pull2(k+1,pl,p2);
       m1=vl[k]; m2=vl[k+1];
       if(fabs(m1)<Eps){
         addptd2(ql, p1);
         continue;
       }
       if(fabs(m2)<Eps){
         continue;
       }
       if((m1>0)&&(m2>0)){
         continue;
       }
       if((m1< 0)&&(m2< 0)){
         continue;
       }
       tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
       tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
       addptd2(ql, tmpd2);
     }
     uval1=uval2;
     eval11=eval21;
     eval12=eval22;
     if(length2(ql)==2){
       puv1[0]=ql[1][0]; puv1[1]=ql[1][1];
       puv2[0]=ql[2][0]; puv2[1]=ql[2][1];
       surffun(ch,puv1[0],puv1[1],p1d3);
       surffun(ch,puv2[0],puv2[1],p2d3);
       v1=Vec[0]; v2=Vec[1]; v3=Vec[2];
       if(fabs(v1)>Eps){
         m1=pta[2]+v3/v1*(p1d3[0]-pta[0])-p1d3[2];
         m2=pta[2]+v3/v1*(p2d3[0]-pta[0])-p2d3[2];
       }else if(fabs(v2)>Eps){
         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[2];
       }
       if(m1*m2>= 0){
         if(((m1>0) && (m2>0)) || ((m1< 0) && (m2< 0))){
           continue;
         }
         if(m1==0){
           ptd3[0]=p1d3[0]; ptd3[1]=p1d3[1]; ptd3[2]=p1d3[2];
           ptuv[0]=puv1[0]; ptuv[1]=puv1[1];
         }else{
           ptd3[0]=p2d3[0]; ptd3[1]=p2d3[1]; ptd3[2]=p2d3[2];
           ptuv[0]=puv2[0]; ptuv[1]=puv2[1];
         }
       }else{
         ptd3[0]=1/(m1-m2)*(-m2*p1d3[0]+m1*p2d3[0]);
         ptd3[1]=1/(m1-m2)*(-m2*p1d3[1]+m1*p2d3[1]);
         ptd3[2]=1/(m1-m2)*(-m2*p1d3[2]+m1*p2d3[2]);
         ptuv[0]=1/(m1-m2)*(-m2*puv1[0]+m1*puv2[0]);
         ptuv[1]=1/(m1-m2)*(-m2*puv1[1]+m1*puv2[1]);
       }
       tmpd3[0]=ptd3[0]-pta[0];
       tmpd3[1]=ptd3[1]-pta[1];
       tmpd3[2]=ptd3[2]-pta[2];
       crossprod(3,tmpd3,Vec, tmp1d3);
       if(norm(3,tmp1d3)<Eps1){
                 tmpd1=zparapt(ptd3)-zparapt(pta)-Eps2;
         if(tmpd1>0){
               copyd(0,2,ptd3,out); out[3]=tmpd1;
           return 1;
         }else{
           if(tmpd1>out[3]){
                 copyd(0,2,ptd3,out); out[3]=tmpd1;
           }
         }
       }
     }
   }
 }
 return 0;
}

int nohiddenpara2(short ch,double par[], double fk[][3], int uveq, double figkL[][3]){
 double eps0=pow(10,-4.0), s, p1[2], p2[2], q[2],tmpd1, tmpd2[2],tmp1,tmp2;
 double paL[DsizeM], fh[DsizeM][2], pta[3], ptap[2], vec[3];
 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 flg;
 int ncusp=floor(Cusppt[0][0]+0.5), seL[DsizeM];
 int i,j,n, nfh, nfk, cspflg,tmpi1[DsizeM];
 vec[0]=sin(THETA)*cos(PHI);
 vec[1]=sin(THETA)*sin(PHI);
 vec[2]=cos(THETA);
 nfk=length3(fk);
 nfh=projpara(fk, fh);
 pull2(1,fh,p1);
 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;
   }
 }
 paL[0]=0;
 tmp2=-1;
 for(j=1;j<=length1(par)-1;j++){
   tmp1=par[j];
   if(fabs(tmp1-tmp2)>Eps){
     add1(paL,tmp1);
     tmp2=tmp1;
   }
 }
 tmp1=paL[length1(paL)];
 tmp2=par[length1(par)];
 if(fabs(tmp1-tmp2)<Eps){
   paL[length1(paL)]=tmp2;
 }else{
   add1(paL,tmp2);
 }
 seL[0]=0;
 for(n=1; n<=length1(paL)-1; n++){
       s=(paL[n]+paL[n+1])/2.0;
   invparapt(s,fh,fk, pta);
   parapt(pta,ptap);
   flg=pthiddenQ(ch,pta, uveq, tmpd4);
   if(flg==0){ //180518
     appendi1(n,seL);
   }
 }
 figL[0][0]=0;figkL[0][0]=0;
 for(i=1; i<=seL[0]; i++){
   sp=paL[seL[i]]; //180616from
   sq=paL[seL[i]+1];
   pointoncurve(sp, fh, Eps, pd2);
   pointoncurve(sq, fh, Eps, qd2);
   if(dist(2,pd2,qd2)>Eps1){
     nearestptpt(pd2,fh,tmpd4);
     dtmp1=tmpd4[2];
     nearestptpt(qd2,fh,tmpd4);
     dtmp2=tmpd4[2];
     partcrv(dtmp1,dtmp2,fh,tmpmd2);
     appendd2(2,1,length2(tmpmd2),tmpmd2,figL);
     tp=invparapt(sp,fh,fk,tmpd3);
     tq=invparapt(sq,fh,fk,tmpd3);
     partcrv3(tp,tq,fk,tmpmd3);
     appendd3(2,1,length3(tmpmd3), tmpmd3, figkL);
   }//180616to
   else{//180705from
     if(fabs(sp-sq)>1){
       appendd2(2,1,length2(fh),fh,figL);
       appendd3(2,1,length3(fk),fk,figkL);
     }
   }//180705to
 }
 tmpi1[0]=0;
 for(i=1; i<=length1(paL)-1; i++){
   if(memberi(i,seL)==0){
     appendi1(i,tmpi1);
   }
 }
 copyi(0,tmpi1[0],tmpi1,seL);
 seL[0]=tmpi1[0];
 Hiddendata[0][0]=0;
 for(i=1; i<=seL[0]; i++){
   dtmp=paL[seL[i]];
   pointoncurve(dtmp,fh,Eps,tmp1d2);
   dtmp=paL[seL[i]+1];
   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);
       partcrv3(tp,tq,fk,tmpmd3);
       appendd3(2,1,length3(tmpmd3),tmpmd3,Hiddendata);
       copyd(0,1,tmp1d2,pd2);
       sp=paL[seL[i]];
       copyd(0,1,tmp2d2,qd2);
       sq=paL[seL[i]+1];
     }
   }
 }
 if(seL[0]>0){
   if(dist(2,pd2,qd2)>Eps1){
     tp=invparapt(sp,fh,fk,tmpd3);
     tq=invparapt(sq,fh,fk,tmpd3);
     partcrv3(tp,tq,fk,tmpmd3);
     appendd3(2,1,length3(tmpmd3),tmpmd3,Hiddendata);
   }
   else{
     appendd3(2,1,length3(fk),fk,Hiddendata);
   }
 }
 return length3(figkL);
}

void borderparadata(short ch,double fkL[][3], double fsL[][3]){
 double ekL[DsizeL][3], fall[DsizeL][2], fbdxy[DsizeL][2];
 double tmpd3[DsizeM][3],  tmpd2[DsizeM][2], p[2], p3[3];
 double figkL[DsizeM][3], du, dv;
 int i, j;
 int din[DsizeS][2],din2[DsizeM][2];
 int ntmp, ntmpd3,ntmpd2;
 double par[DsizeM];
 du=(Urng[1]-Urng[0])/Mdv;
 dv=(Vrng[1]-Vrng[0])/Ndv;
 fsL[0][0]=0;
 ekL[0][0]=0;;
 if(DrawE==1){
   tmpd2[0][0]=0;;
   for(j=0; j<=Ndv; j++){
     add2(tmpd2, Urng[1], Vrng[0]+j*dv);
   }
   surfcurve(ch,tmpd2,tmpd3);
   appendd3(2,1,length3(tmpd3),tmpd3,ekL);
 }
 if(DrawN==1){
   tmpd2[0][0]=0;
   for(j=0; j<=Mdv; j++){
     add2(tmpd2, Urng[0]+j*du, Vrng[1]);
   }
   surfcurve(ch,tmpd2,tmpd3);
   appendd3(2,1,length3(tmpd3),tmpd3,ekL);
 }
 if(DrawW==1){
   tmpd2[0][0]=0;
   for(j=0; j<=Ndv; j++){
     add2(tmpd2, Urng[0], Vrng[0]+j*dv);
   }
   surfcurve(ch,tmpd2,tmpd3);
   appendd3(2,1,length3(tmpd3),tmpd3,ekL);
 }
 if(DrawS==1){
   tmpd2[0][0]=0;
   for(j=0; j<=Mdv; j++){
     add2(tmpd2, Urng[0]+j*du, Vrng[0]);
   }
   surfcurve(ch,tmpd2,tmpd3);
   appendd3(2,1,length3(tmpd3),tmpd3,ekL);
 }
 if(length3(ekL)>0){
   appendd3(2,1,length3(ekL),ekL,fkL);
 }
 projpara(fkL,fall);
 makexybdy(ch,fbdxy);
 dataindexd2(2,fall,din);
 Borderpt[0][0]=0;
 Otherpartition[0]=0;
 for(i=1; i<=length2i(din)-1; i++){
   Otherpartition[i]=Inf; Otherpartition[0]++;
 }
 Borderhiddendata[0][0]=0;
 dataindexd3(2,fkL,din2);
 for(i=1; i<= length2i(din2); i++){
       tmpd3[0][0]=0;
   appendd3(0,din2[i][0],din2[i][1],fkL,tmpd3);
   projpara(tmpd3,tmpd2);
       partitionseg(tmpd2,fall,0,par);
/*
       if(length1(par)>2){
     ntmp=length1(par)-1;
     tmpd2[0][0]=0;
     for(j=2; j<=ntmp; j++){
       pull2(j,Partitionpt,p);
       addptd2(tmpd2, p);
     }
     appendd2(0, 1,length2(tmpd2),tmpd2,Borderpt);
   }
*/
       tmpd3[0][0]=0;
   appendd3(0,din2[i][0],din2[i][1],fkL,tmpd3);
//    for(j=din2[i][0]; j<=din2[i][1]; j++){
//      pull3(j,fkL,p3);
//      addptd3(tmpd3,p3);
//    }
       nohiddenpara2(ch,par,tmpd3,1, figkL);
       if(length3(figkL)>0){
     appendd3(2,1,length3(figkL),figkL, fsL);
   }
       ntmp=length3(Hiddendata);
   if(ntmp>0){
     appendd3(2,1,ntmp,Hiddendata,Borderhiddendata);
   }
 }
}

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(1,ptL);
   pull3(1,pd,p);
   for(k=2; k<=npd-1; k++){
     pull3(k,pd,tmp2d);
     if(dist(3,p,tmp2d)>eps1){
       nptL=appendi1(k,ptL);
       pull3(k,pd,p);
     }
   }
   pull3(npd-1,pd,p);
   pull3(npd,pd,tmp2d);
   if(dist(3,p,tmp2d)>eps){  // eps -> eps1 ?
     nptL=appendi1(npd,ptL);
   }
   if(nout>0){
     nout=appendi1(Infint,out);
   }
   for(i=1; i<=nptL; i++){
     nout=appendi1(ptL[i],out);
   }
 }
 out[0]=nout;
 return nout;
}

void sfbdparadata(short ch,double outd3[][3]){
 double pdatad3[DsizeL][3];
 double pts[DsizeL][2],out3md[DsizeL][2],eps;
 double tmpmd[DsizeL][2],tmp1md[DsizeL][2],tmp2md[DsizeL][2],tmp3md[DsizeL][2];
 double tmpd[2],tmp1d[2],tmpd3[3],tmp1d3[3],tmp2d3[3],tmp2md3[DsizeL][3];
 int din[DsizeS][2],nlist[DsizeL];
 int ii,jj,kk,n,nall,flg;
 envelopedata(ch, out3md); //180517from
 tmp3md[0][0]=0;
 pts[0][0]=0;
 dataindexd2(2,out3md,din);
 for(jj=1;jj<=length2i(din);jj++){
       tmp1md[0][0]=0;
   tmp2md3[0][0]=0;
   for(kk=din[jj][0];kk<=din[jj][1];kk++){
     pull2(kk,out3md,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,tmp3md);
     }
   }
 }
 out3md[0][0]=0;
 appendd2(0,1,length2(tmp3md),tmp3md,out3md);
 tmp3md[0][0]=0;
 appendd2(0,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++){
   pull2(ii,pts,tmpd);
   surffun(ch,tmpd[0],tmpd[1],tmp1d3);
   dataindexd2(2,out3md,din);
   for(jj=1;jj<=length2i(din);jj++){
     tmp2md[0][0]=0;
     appendd2(0,din[jj][0],din[jj][1],out3md,tmp2md);
     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){
       prependd2(ii,ii,pts,tmp2md);
       replacelistd2(2,jj,out3md,tmp2md);
       dataindexd2(2,out3md,din);//180517
       continue;
     }
     pull2(length2(tmp2md),tmp2md,tmpd);
     surffun(ch,tmpd[0],tmpd[1],tmpd3);
     if(dist(3,tmp1d3,tmpd3)<eps+Eps){
       appendd2(0,ii,ii,pts,tmp2md);
       replacelistd2(2,jj,out3md,tmp2md);
       dataindexd2(2,out3md,din);//180517
       continue;
     }
   }
 }
 nall=cuspsplitpara(ch,out3md, pdatad3);
 n=length2(Cuspsplitpt);
 Cusppt[0][0]=0;
 nall=appendd2(0, 1,n,Cuspsplitpt,Cusppt);
 borderparadata(ch,pdatad3, outd3);
 push3(Inf,3,0,outd3);
 nall=appendd3(0,1,length3(Borderhiddendata),Borderhiddendata,outd3);
 push3(Inf,3,0,outd3);
}

void groupnearpt3(double qlist[][3],double plist[][3]){
// 180428
 int ii,jj,flg,din[DsizeS][2];
 double p1[3],p2[3],gl[DsizeS][3];
 while(length3(qlist)>0){
   gl[0][0]=0;
   appendd3(2,1,1,qlist,gl);
   removed3(1,qlist);
   for(ii=1;ii<=length3(qlist);ii++){
     flg=0;
     pull3(ii,qlist,p1);
     for(jj=1;jj<=length3(gl);jj++){
       pull3(jj,gl,p2);
       if(dist(3,p1,p2)<Eps1){
         flg=1;
         break;
       }
     }
     if(flg==1){
       addptd3(gl,p1);
       removed3(ii,qlist);
     }
   }
   appendd3(2,1,length3(gl),gl,plist);
 }
}

double funmeet(short ch,double u, double v,double pa[3], double vec[3]){
 double pt[3], out;
 surffun(ch,u, v, pt);
 if((fabs(vec[1])>Eps)||(fabs(vec[0])>Eps)){
   out=(vec[1]*(pt[0]-pa[0])-vec[0]*(pt[1]-pa[1]));
 }else{
   out=pt[0];
 }
 return out;
}

void meetpoints(short ch,double pta[3], double ptb[3], int uveq,double out[][3]){
 int i, j, k,flg,din[DsizeS][2];
 double vec[3],v1,v2,v3;
 double du=(Urng[1]-Urng[0])/Mdv, dv=(Vrng[1]-Vrng[0])/Ndv;
 double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
 double pl[5][2], vl[5], ql[11][2],p1[2],p2[2],m1,m2;
 double ptuv[2], puv1[2], puv2[2],p1d3[3], p2d3[3], ptd3[3];
 double tmpd1, tmpd2[2],tmpd3[3],tmp1d3[3],tmp1md3[DsizeS][3];
 vec[0]=ptb[0]-pta[0];  vec[1]=ptb[1]-pta[1];  vec[2]=ptb[2]-pta[2];
 out[0][0]=0;
 if(norm(3,vec)<Eps){
   return;
 }
 for(j=0; j<=Ndv-1; j++){
   vval1=Vrng[0]+j*dv;
   vval2=Vrng[0]+(j+1)*dv;
   uval1=Urng[0];
   eval11=funmeet(ch,uval1,vval1,pta,vec);
   eval12=funmeet(ch,uval1,vval2,pta,vec);
   for(i=0; i<=Mdv-1; i++){
     uval2=Urng[0]+(i+1)*du;
     eval21=funmeet(ch,uval2,vval1,pta,vec);
     eval22=funmeet(ch,uval2,vval2,pta,vec);
         pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
     pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
     pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
     pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
         pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
     ql[0][0]=0;
     for(k=0; k<=3; k++){
       pull2(k,pl,p1);
       pull2(k+1,pl,p2);
       m1=vl[k]; m2=vl[k+1];
       if(fabs(m1)<Eps){
         addptd2(ql, p1);
         continue;
       }
       if(fabs(m2)<Eps){
         continue;
       }
       if((m1>0)&&(m2>0)){
         continue;
       }
       if((m1< 0)&&(m2< 0)){
         continue;
       }
       tmpd2[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
       tmpd2[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
       addptd2(ql, tmpd2);
     }
     uval1=uval2;
     eval11=eval21;
     eval12=eval22;
     if(length2(ql)==2){
       puv1[0]=ql[1][0]; puv1[1]=ql[1][1];
       puv2[0]=ql[2][0]; puv2[1]=ql[2][1];
       surffun(ch,puv1[0],puv1[1],p1d3);
       surffun(ch,puv2[0],puv2[1],p2d3);
       v1=vec[0]; v2=vec[1]; v3=vec[2];
       if(fabs(v1)>Eps){
         m1=pta[2]+v3/v1*(p1d3[0]-pta[0])-p1d3[2];
         m2=pta[2]+v3/v1*(p2d3[0]-pta[0])-p2d3[2];
       }else if(fabs(v2)>Eps){
         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[2];
       }
       if(m1*m2>= 0){
         if(((m1>0) && (m2>0)) || ((m1< 0) && (m2< 0))){
           continue;
         }
         if(m1==0){
           ptd3[0]=p1d3[0]; ptd3[1]=p1d3[1]; ptd3[2]=p1d3[2];
           ptuv[0]=puv1[0]; ptuv[1]=puv1[1];
         }else{
           ptd3[0]=p2d3[0]; ptd3[1]=p2d3[1]; ptd3[2]=p2d3[2];
           ptuv[0]=puv2[0]; ptuv[1]=puv2[1];
         }
       }else{
         ptd3[0]=1/(m1-m2)*(-m2*p1d3[0]+m1*p2d3[0]);
         ptd3[1]=1/(m1-m2)*(-m2*p1d3[1]+m1*p2d3[1]);
         ptd3[2]=1/(m1-m2)*(-m2*p1d3[2]+m1*p2d3[2]);
         ptuv[0]=1/(m1-m2)*(-m2*puv1[0]+m1*puv2[0]);
         ptuv[1]=1/(m1-m2)*(-m2*puv1[1]+m1*puv2[1]);
       }
       tmpd3[0]=ptd3[0]-pta[0];
       tmpd3[1]=ptd3[1]-pta[1];
       tmpd3[2]=ptd3[2]-pta[2];
       crossprod(3,tmpd3,vec, tmp1d3);
       if(norm(3,tmp1d3)<Eps1){
         tmpd1=dotprod(3,tmpd3,vec);
         if((tmpd1>-Eps)&&(norm(3,tmpd3)<norm(3,vec)+Eps)){
           if(length3(out)==0){
             addptd3(out,ptd3);
           }else{
             pull3(length3(out),out,p1d3);
             if(dist(3,p1d3,ptd3)>Eps){
               addptd3(out,ptd3);
             }
           }
         }
       }
     }
   }
 }
 tmp1md3[0][0]=0;
 groupnearpt3(out,tmp1md3);
 dataindexd3(2,tmp1md3,din);
 out[0][0]=0;
 for(i=1;i<=din[0][0];i++){
   p1d3[0]=0; p1d3[1]=0; p1d3[2]=0;
   tmpd1=din[i][1]-din[i][0]+1;
   for(j=din[i][0];j<=din[i][1];j++){
     p1d3[0]=p1d3[0]+tmp1md3[j][0]/tmpd1;
     p1d3[1]=p1d3[1]+tmp1md3[j][1]/tmpd1;
     p1d3[2]=p1d3[2]+tmp1md3[j][2]/tmpd1;
   }
   addptd3(out,p1d3);
 }
}

void crvsfparadata(short ch, double fkL[][3], double fbdkL[][3], int sepflg, double out[][3]){
 double fbdy[DsizeL][2], fk[DsizeM][3],outh[DsizeL][3];
//  double fkp[DsizeM][3], dt;
 double fh[DsizeM][2], parL[DsizeM], tmpmd3[DsizeM][3];
 double ptL[DsizeM][3], tmpd2[2], tmpd3[3], tmp1d1, tmp2d1;
 double po[2]={0,0}, epsd2[2]={Eps1,1}, pa[3], pb[3];
 int nbor=length3(Borderhiddendata);
 int ncshidden,din[DsizeM][2],din2[DsizeS][2];
 int nn, i, j, k, n;
 projpara(fbdkL,fbdy);
 out[0][0]=0;
 outh[0][0]=0;
 dataindexd3(2,fkL,din);
 for(nn=1; nn<=length2i(din); nn++){
   fk[0][0]=0;
   appendd3(0,din[nn][0],din[nn][1],fkL,fk);
   projpara(fk,fh);
   partitionseg(fh,fbdy,0,parL);
   if(sepflg==0){ //180522
     for(i=1; i<=length3(fk)-1; i++){
       pull3(i,fk,pa);
       pull3(i+1,fk,pb);
       meetpoints(ch,pa,pb,1,ptL);
       for(j=1; j<=length3(ptL); j++){
             pull3(j,ptL,tmpd3);
         parapt(tmpd3,tmpd2);
         tmp1d1=paramoncurve(tmpd2,i,fh);
         tmp2d1=Inf;
         for(k=1; k<=length1(parL); k++){
            tmp2d1=fmin(tmp2d1,fabs(parL[k]-tmp1d1));
         }
         tmpd3[0]=pb[0]-pa[0]; tmpd3[1]=pb[1]-pa[1];
         tmpd3[2]=pb[2]-pa[2];
         if(tmp2d1*dist(3,tmpd3,po)>Eps1){
           appendd1(tmp1d1,parL);
         }
       }
       simplesort(parL);
     }
   }
   nohiddenpara2(ch,parL,fk,1, tmpmd3);
   appendd3(2,1,length3(tmpmd3),tmpmd3,out);
   appendd3(2,1, length3(Hiddendata),Hiddendata,outh);
 }
 push3(Inf,3,0,out); //17.06.17
 appendd3(0,1,length3(outh),outh,out);
}

void crv3onsfparadata(short ch,double fk[][3], double fbdyd3[][3], double out[][3]){
 // 180609 debugged(renewed)
 double fbdy[DsizeL][2],fh[DsizeL][2],fks[DsizeL][3],fhs[DsizeL][2],par[DsizeM];
 double tmpmd3[DsizeL][3],outh[DsizeL][3];
 int i,j,din[DsizeS][2];
 projpara(fbdyd3,fbdy);
 projpara(fk,fh);
 out[0][0]=0;
 outh[0][0]=0;
 dataindexd2(2,fh,din);
 for(i=1;i<=din[0][0];i++){
   fhs[0][0]=0; fks[0][0]=0;
       appendd2(0,din[i][0],din[i][1],fh,fhs);
   appendd3(0,din[i][0],din[i][1],fk,fks);
   tmpmd3[0][0]=0;
   if(length2(fhs)>1){
     partitionseg(fhs,fbdy,0, par);
     nohiddenpara2(ch,par,fks,1,tmpmd3);
     appendd3(2,1,length3(tmpmd3),tmpmd3,out);
     appendd3(2,1,length3(Hiddendata),Hiddendata,outh);
   }else{
     appendd3(2,1,1,fks,out);
   }
 }
 connectseg3(out, Eps1,tmpmd3);
 out[0][0]=0;
 appendd3(0,1,length3(tmpmd3),tmpmd3,out);
 push3(Inf,3,0,out);  //181025
 connectseg3(outh, Eps1,tmpmd3);
 appendd3(0,1,length3(tmpmd3),tmpmd3,out);//181025
}

void crv2onsfparadata(short ch,double fh[][2], double fbdyd3[][3], double out[][3]){
 double fk[DsizeL][3];
 surfcurve(ch,fh,fk);
 crv3onsfparadata(ch,fk,fbdyd3,out);
}

void wireparadata(short ch,double bdyk[][3], double udata[], double vdata[],const char *fname,const char *fnameh){
 double crv2d[DsizeL][2],out[DsizeL][3],out1[DsizeL][3],out2[DsizeL][3];
 double du=(Urng[1]-Urng[0])/Mdv;
 double dv=(Vrng[1]-Vrng[0])/Ndv;
 int i,j,din[DsizeS][2];
 short allflg=1;
 if(strlen(fnameh)>0){
   allflg=0;
 }
 char var0[]="wire3d";
 char varh0[]="wireh3d";
 char var[20];
 char varh[20];
 char varnow[40];
 char varhnow[40];
 varnow[0]='\0';
 varhnow[0]='\0';
 sprintf(var,"%s%d",var0,ch);
 sprintf(varh,"%s%d",varh0,ch);
 char dirfname[256];
 char dirfnameh[256];
 char varname[256];
 char varnameh[256];
 dirfname[0]='\0';
 dirfnameh[0]='\0';
 varname[0]='\0';
 varnameh[0]='\0';
 sprintf(varname,"%s%d",var,ch);
 sprintf(varnameh,"%s%d",varh,ch);
 sprintf(dirfname,"%s%s",Dirname,fname);
 sprintf(dirfnameh,"%s%s",Dirname,fnameh);
 if((length1(udata)==0)&&(allflg==0)){
   out[0][0]=0;
   output3(1,"w",varname,dirfname,out);
   output3(1,"w",varnameh,dirfnameh,out);
 }
 for(i=1;i<=length1(udata);i++){
   crv2d[0][0]=0;
   for(j=0; j<=Ndv; j++){
     add2(crv2d, udata[i],Vrng[0]+j*dv);
   }
   if(allflg==1){out[0][0]=0;}
       crv2onsfparadata(ch,crv2d,bdyk,out);
   if(allflg==0){
     dataindexd3(3,out,din);
     out1[0][0]=0; out2[0][0]=0;
     appendd3(0,din[1][0],din[1][1],out,out1);
     appendd3(0,din[2][0],din[2][1],out,out2);
     if(i==1){
       output3(1,"w",varname,dirfname,out1);
       output3(1,"w",varnameh,dirfnameh,out2);
     }else{
       output3(0,"a",varname,dirfname,out1);
       output3(0,"a",varnameh,dirfnameh,out2);
     }
   }else{
     sprintf(varnow,"%s%s%d",var,"u",i);
     sprintf(varhnow,"%s%s%d",varh,"u",i);
     output3h("a",varnow, varhnow,fname,out);
   }
 }
 for(j=1;j<=length1(vdata);j++){
   crv2d[0][0]=0;
   for(i=0; i<=Mdv; i++){
     add2(crv2d, Urng[0]+i*du,vdata[j]);
   }
   if(allflg==1){out[0][0]=0;}
   crv2onsfparadata(ch,crv2d,bdyk,out);
   if(allflg==0){
     dataindexd3(3,out,din);
     out1[0][0]=0; out2[0][0]=0;
     appendd3(0,din[1][0],din[1][1],out,out1);
     appendd3(0,din[2][0],din[2][1],out,out2);
     if(i<length1(vdata)){
       output3(0,"a",varname,dirfname,out1);
       output3(0,"a",varnameh,dirfnameh,out2);
     }else{
       output3(0,"a",varname,dirfname,out1);
       outputend(dirfname);
       output3(0,"a",varnameh,dirfnameh,out2);
       outputend(dirfnameh);
     }
   }else{
     sprintf(varnow,"%s%s%d",var,"v",j);
     sprintf(varhnow,"%s%s%d",varh,"v",j);
         output3h("a",varnow, varhnow,fname,out);
   }
 }
 if(allflg==0){
   outputend(dirfname);
   outputend(dirfnameh);
 }
}

void intersectcrvsf(const char *wa, short chfd,double crv[][3],const char *fname){
 double pa[3],pb[3],out[DsizeM][3],ptL[DsizeM][3];
 int i;
 short chcut;
 char chc[10];
 sprintf(chc,"%d",chfd);
 char var[]="intercrvsf";
 char dirfname[256];
 dirfname[0] = '\0';
 sprintf(dirfname,"%s%s",Dirname,fname);
 char varname[256];
 varname[0] = '\0';
 sprintf(varname,"%s%s",var,chc);
 ptL[0][0]=0;
 for(i=1;i<length3(crv);i++){
   pull3(i,crv,pa);
   pull3(i+1,crv,pb);
   meetpoints(chfd,pa,pb,1,out);
   appendd3(0,1,length3(out),out,ptL);
 }
 output3(1,wa,varname,dirfname,ptL);
}

void sfcutparadata(short chfd, short ncut, double fbdy3[][3],const char *fname,const char *fnameh){
 double Vec[3]={sin(THETA)*cos(PHI),sin(THETA)*sin(PHI),cos(THETA)};
 double v1=Vec[0],v2=Vec[1],v3=Vec[2];
 double out[DsizeM][2],out2[DsizeM][2],outd3[DsizeL][3];
 double ekL[DsizeL][3],out1d3[DsizeL][3],out2d3[DsizeL][3];
 double uval1,uval2,vval1,vval2,eval11,eval12,eval21,eval22;
 double pl[5][2], vl[5], ql[11][2],diff1,diff2;
 double tmp1md[DsizeM][2],tmp2md[DsizeM][2],tmpd[2],tmp1d[2],tmp2d[2];
 double tmp2md3[DsizeM][3],tmpd3[3],tmp1d3[3];
 double du=(Urng[1]-Urng[0])/Mdv;
 double dv=(Vrng[1]-Vrng[0])/Ndv;
 double p1[2],p2[2],m1,m2,q11,q12,q21,q22,a1,a3,b1,b3;
 int i,j,k,din[DsizeS][2],ns,ne,rmv[DsizeS];
 short chcut, allflg=0;
 if(strlen(fnameh)==0){
   allflg=1;
 }
 char var0[]="sfcut3d";
 char varh0[]="sfcuth3d";
 char var[20];
 char varh[20];
 sprintf(var,"%s%d",var0,chfd);
 sprintf(varh,"%s%d",varh0,chfd);
 char varnow[40];
 char varhnow[40];
 char dirfname[256];
 char dirfnameh[256];
 char varname[256];
 char varnameh[256];
 varnow[0]='\0';
 varhnow[0]='\0';
 dirfname[0]='\0';
 dirfnameh[0]='\0';
 varname[0]='\0';
 varnameh[0]='\0';
 sprintf(varname,"%s%d",var,chfd);
 sprintf(varnameh,"%s%d",varh,chfd);
 sprintf(dirfname,"%s%s",Dirname,fname);
 sprintf(dirfnameh,"%s%s",Dirname,fnameh);
 for(chcut=1;chcut<=ncut;chcut++){
   out[0][0]=0; out2[0][0]=0;
   for(j=0; j<=Ndv-1; j++){
     vval1=Vrng[0]+j*dv;
     vval2=Vrng[0]+(j+1)*dv;
     uval1=Urng[0];
     eval11=cutfun(chfd, chcut,uval1, vval1);
     eval12=cutfun(chfd, chcut,uval1, vval2);
     for(i=0; i<=Mdv-1; i++){
       uval2=Urng[0]+(i+1)*du;
       eval21=cutfun(chfd, chcut,uval2, vval1);
       eval22=cutfun(chfd, chcut,uval2, vval2);
           pl[0][0]=uval1; pl[0][1]=vval1; vl[0]=eval11;
       pl[1][0]=uval2; pl[1][1]=vval1; vl[1]=eval21;
       pl[2][0]=uval2; pl[2][1]=vval2; vl[2]=eval22;
       pl[3][0]=uval1; pl[3][1]=vval2; vl[3]=eval12;
           pl[4][0]=uval1; pl[4][1]=vval1; vl[4]=eval11;
       ql[0][0]=0;
       for(k=0; k<=3; k++){
         pull2(k,pl,p1);
         pull2(k+1,pl,p2);
         m1=vl[k]; m2=vl[k+1];
         if(fabs(m1)<Eps){
           addptd2(ql, p1);
           continue;
         }
         if(fabs(m2)<Eps){
           continue;
         }
         if((m1>0)&&(m2>0)){
           continue;
         }
         if((m1< 0)&&(m2< 0)){
           continue;
         }
         tmpd[0]=1/(m1-m2)*(-m2*p1[0]+m1*p2[0]);
         tmpd[1]=1/(m1-m2)*(-m2*p1[1]+m1*p2[1]);
         addptd2(ql, tmpd);
       }
       uval1=uval2;
       eval11=eval21;
       eval12=eval22;
       if(length2(ql)==2){
         q11=ql[1][0];q12=ql[1][1];q21=ql[2][0];q22=ql[2][1];
         a1=uval1;a3=uval2;b1=vval1;b3=vval2;
         if(((q11==a1)&&(q21==a1))||((q11==a3)&&(q21==a3))){
           if(((q21==b1)&&(q22==b1))||((q21==b3)&&(q22==b3))){
             appendd2(2,1,2,ql,out2);
           }else{
             appendd2(2,1,2,ql,out);
           }
           continue;
         }
         if(((q12==b1)&&(q22==b1))||((q12==b3)&&(q22==b3))){
           if(((q11==a1)&&(q21==a1))||((q11==a3)&&(q21==a3))){
             appendd2(2,1,2,ql,out2);
           }else{
             appendd2(2,1,2,ql,out);
           }
           continue;
         }
         appendd2(2,1,2,ql,out);
       }
     }
   }
   dataindexd2(2,out2,din);
   while(length2i(din)>0){
     ns=din[1][0];ne=din[1][1];
     tmp1md[0][0]=0;
     appendd2(0,ns,ne,out2,tmp1md);
     appendd2(2,ns,ne,out2,out);
     removelistd2(2,1,out2);
     dataindexd2(2,out2,din);
     rmv[0]=0;
     for(j=1;j<=length2i(din);j++){
       ns=din[j][0];ne=din[j][1];
       tmp2md[0][0]=0;
       appendd2(0,ns,ne,out2,tmp2md);
       pull2(1,tmp2md,tmp2d);pull2(1,tmp1md,tmp1d);
       diff1=dist(2,tmp1d,tmp2d);
       pull2(2,tmp2md,tmp2d);pull2(2,tmp1md,tmp1d);
       diff1=diff1+dist(2,tmp1d,tmp2d);
       pull2(1,tmp2md,tmp2d);pull2(2,tmp1md,tmp1d);
       diff2=dist(2,tmp1d,tmp2d);
       pull2(2,tmp2md,tmp2d);pull2(1,tmp1md,tmp1d);
       diff2=diff2+dist(2,tmp1d,tmp2d);
       if((diff1<Eps)||(diff2<Eps)){
         appendi1(j,rmv);
         continue;
       }
     }
     for(j=1;j<=rmv[0];j++){
       removelistd2(2,rmv[j],out2);
     }
     dataindexd2(2,out2,din);
   }
   connectseg(out, Eps1, out2);
   dataindexd2(2,out2,din);
   outd3[0][0]=0;
   for(i=1;i<=length2i(din);i++){
     ns=din[i][0];ne=din[i][1];
     tmp1md[0][0]=0;
     appendd2(0,ns,ne,out2,tmp1md);
     tmp2md3[0][0]=0;
     for(j=1;j<=length2(tmp1md);j++){
       pull2(j,tmp1md,tmpd);
       surffun(chfd,tmpd[0],tmpd[1],tmpd3);
       if(length3(tmp2md3)>0){
         pull3(length3(tmp2md3),tmp2md3,tmp1d3);
         if(dist(3,tmpd3,tmp1d3)<Eps1){
           continue;
         }
       }
       addptd3(tmp2md3,tmpd3);
     }
     appendd3(2,1,length3(tmp2md3),tmp2md3,outd3);
       }
   crv3onsfparadata(chfd,outd3,fbdy3,ekL);
   if(allflg==0){
     dataindexd3(3,ekL,din);
     out1d3[0][0]=0; out2d3[0][0]=0;
     appendd3(0,din[1][0],din[1][1],ekL,out1d3);
     appendd3(0,din[2][0],din[2][1],ekL,out2d3);
     if(chcut==1){
       output3(1,"w",varname,dirfname,out1d3);
       output3(1,"w",varnameh,dirfnameh,out2d3);
     }else{
       output3(0,"a",varname,dirfname,out1d3);
       output3(0,"a",varnameh,dirfnameh,out2d3);
     }
   }else{
     sprintf(varnow,"%s%d",var,chcut);
     sprintf(varhnow,"%s%d",varh,chcut);
     output3h("a",varnow, varhnow,fname,ekL);
   }
 }
 if(allflg==0){
   outputend(dirfname);
   outputend(dirfnameh);
 }
}

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 readoutdata3(const char *fname, const char *var, double data[][3]){
 double x,y,z;
 float xx;
 char dstrorg[256];
 dstrorg[0]='\0';
 char dstr[256];
 dstr[0]='\0';
 char str[10];
 str[0]='\0';
 char strtmp[30];
 char tmp[10];
 tmp[0]='\0';
 int linectr=0, start=0,ii,jj,nn,nctr;
 FILE *fp;
 fp=fopen(fname,"r");
 if(fp==NULL){
   printf("file not found\n");
   return;
 }
 data[0][0]=0;
 nn=strlen(var);
 while( !feof(fp)){
   fgets(dstr,250,fp);
   linectr++;
   jj=strlen(dstr);
   if(jj>1){
     dstr[jj-3]='\0';
   }else{
     dstr[0]='\0';
   }
   if(strncmp(dstr,var,nn)==0){
     start=linectr;
   }else{
     if(start==0){
       continue;
     }
   }
   if(linectr==start){
     continue;
   }
   if(strncmp(dstr,"sta",3)==0){
     if(linectr>start+1){
       add3(data,Inf,2,0);
     }
     continue;
   }
   if(strncmp(dstr,"end",3)==0){
     continue;
   }
   if(strncmp(dstr,"[",1)!=0){
     break;
   }
   str[0]='\0';
   nctr=0;
   for(jj=2;jj<250;jj++){
     tmp[0]='\0'; sprintf(tmp,"%c",dstr[jj]);
     if(strncmp(tmp,"/",1)==0){
       break;
     }
     if(strncmp(tmp," ",1)==0){
       continue;
     }
     if(strncmp(tmp,",",1)==0){
       nctr++;
       if(nctr==1){
         x=atof(str); str[0]='\0';
         continue;
       }
       if(nctr==2){
         y=atof(str);str[0]='\0';
         continue;
       }
     }
     if(strncmp(tmp,"]",1)==0){
       z=atof(str);str[0]='\0';
       add3(data,x,y,z);
       nctr=0;
       jj++;
       tmp[0]='\0'; sprintf(tmp,"%c",dstr[jj]);
       if(strncmp(tmp,",",1)==0){
         jj++;
       }
       continue;
     }
     strtmp[0]='\0';
     sprintf(strtmp,"%s%s",str,tmp);
     str[0]='\0';
     for(ii=0;ii<20;ii++){
       str[ii]=strtmp[ii];
       if(strtmp[ii]=='\0'){
         break;
       }
     }
   }
 }
 fclose(fp);
}

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.