Allink: src/ElPoly/ElPolyRepr.cpp Source File

00001 #include "ElPoly.h"
00002 int ElPoly::ProjectionF(int NBin,int Coord){
00003   if(Coord > 4 || Coord <0) return 1;
00004   int NType = 5;
00005   double *Plot = (double *)calloc(NBin*NBin*NType,sizeof(double));
00006   double InvNBin = 1./(double)NBin;
00007   double RefPos[3] = {0.,0.,0.};
00008   for(int d=0;d<3;d++){
00009     RefPos[d] = Nano->Pos[d]-.5*pEdge(d);
00010   }
00011   if(Coord == 3){
00012     RefPos[0]=pCm(0);RefPos[1]=pCm(1);RefPos[2]=pCm(2);
00013   }
00014   SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
00015   for(int f=NFile[0];f<NFile[1];f++){
00016     Processing(f);
00017     OpenRisk(cFile[f],BF_PART);
00018     ShiftSys(SHIFT_CM);
00019     int NPlot = 0;
00020     //---Projects-against-one-coordinate--
00021     if(Coord < 3){
00022       int coord1 = (Coord+1)%3;
00023       int coord2 = (Coord+2)%3;
00024       for(int p=0;p<pNPart();p++){
00025    double x = pPos(p,coord1) - RefPos[coord1];
00026    x -= floor(x*pInvEdge(coord1))*pEdge(coord1);
00027    double y = pPos(p,coord2) - RefPos[coord2];
00028    y -= floor(y*pInvEdge(coord2))*pEdge(coord2);
00029    int v = (int)(NBin*x*pInvEdge(coord1));
00030    if( v < 0 || v >= NBin) continue;
00031    int vv = (int)(NBin*y*pInvEdge(coord2));
00032    if( vv < 0 || vv >= NBin) continue;
00033    int t = pType(p);
00034    if( t < 0 || t > 3) continue;
00035    if( CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED) )
00036      Plot[(v*NBin+vv)*NType+3] += 1.;
00037    Plot[(v*NBin+vv)*NType+t] += 1.;
00038    if(p<pNPart()-1)
00039      if(pType(p+1) == 1 && pType(p) == 0)
00040        Plot[(v*NBin+vv)*NType+4] += 1.;     
00041       }
00042     }
00043     //---Projects-against-the-radial-coordinate--
00044     else if(Coord == 3){
00045       SetEdge(.5*MAX((pEdge(CLat1)),(pEdge(CLat2))),3);
00046       for(int p=0;p<pNPart();p++){
00047    double x = pPos(p,CLat1) - RefPos[CLat1];
00048    x -= floor(x*pInvEdge(CLat1))*pEdge(CLat1);
00049    double y = pPos(p,CLat2) - RefPos[CLat2];
00050    y -= floor(y*pInvEdge(CLat2))*pEdge(CLat2);
00051    double z = pPos(p,CNorm) - RefPos[CNorm];
00052    z -= floor(z*pInvEdge(CNorm))*pEdge(CNorm);
00053    double r = sqrt(SQR(x)+SQR(y));
00054    int v = (int)(NBin*r*pInvEdge(3));
00055    if( v < 0 || v >= NBin) continue;
00056    int vv = (int)(NBin*pPos(p,CNorm)/pEdge(CNorm));
00057    if( vv < 0 || vv >= NBin) continue;
00058    int t = pType(p);
00059    if( t < 0 || t > 3) continue;
00060    if( CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED) )
00061      Plot[(v*NBin+vv)*NType+3] += 1.;
00062    Plot[(v*NBin+vv)*NType+t] += 1.;
00063    if(p<pNPart()-1)
00064      if(pType(p+1) == 1 && pType(p) == 0)
00065        Plot[(v*NBin+vv)*NType+4] += 1.;     
00066       }
00067     }
00068   }
00069   printf("\n");
00070   //-----writes-the-output-file-------------------
00071   FILE *FileToWrite = NULL;
00072   FileToWrite = fopen("Projection.xyd","w");
00073   fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NBin,ChooseDraw(EL_QUAD1));
00074   int coord1 = (Coord+1)%3;
00075   int coord2 = (Coord+2)%3;
00076   if(Coord == 3){
00077     coord1 = 3;
00078     coord2 = CNorm;
00079   }
00080   double Max[NType];
00081   for(int t=0;t<NType;t++){
00082     Max[t] = Plot[t];
00083     for(int v=0;v<NBin;v++)
00084       for(int vv=0;vv<NBin;vv++)
00085    if(Max[t] < Plot[(v*NBin+vv)*NType+t]) Max[t] = Plot[(v*NBin+vv)*NType+t];
00086     Max[t] = Max[t] <= 0. ? 1. : Max[t];
00087   }
00088   //for(int t=0;t<NType-1;t++){
00089   for(int t=0;t<1;t++){
00090     for(int v=0;v<NBin;v++){
00091       for(int vv=0;vv<NBin;vv++){
00092    int p = (v*NBin+vv)*NType+t;
00093    int c = 0;
00094    if(Plot[p] < .1) continue;
00095    double x = (v)*InvNBin*pEdge(CLat1);
00096    double y = (vv)*InvNBin*pEdge(CLat2);
00097    double dens = Plot[p]/Max[t]*5.+.5*pEdge(CNorm);
00098    double NanoAdded = 0.;//Plot[p]/Max[t]+Plot[((v*NBin+vv)*NType+3]/Max[3];
00099    double Phob = t == 0 ? Plot[(v*NBin+vv)*NType+0]/Max[0] : 0.;
00100    double Phil = t == 1 ? Plot[(v*NBin+vv)*NType+1]/Max[1] : 0.;
00101    fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf)}\n",p,c,t,x,y,dens,NanoAdded,Phob,Phil);
00102       }
00103     } 
00104   }
00105   free(Plot);
00106   fclose(FileToWrite);
00107   return 0;
00108 }
00109 // int ElPoly::LateralDistr(int NBin){
00110 // }
00111 int ElPoly::Surface(int NBin,int Threshold){
00112   double *PlotAv = (double *)calloc(NBin*NBin,sizeof(double));
00113   int NPoint=0;
00114   double AreaAv = 0.,SurfAv = 0.;
00115   double Max=0.;
00116   printf("Edge %lf Threshold %d\n",pEdge(3),Threshold);
00117   for(int f=NFile[0];f<NFile[1];f++){
00118     Processing(f);
00119     if(OpenRisk(cFile[f],BF_PART))return 1;
00120     double Surf=0.,Area =0.;
00121     double *Plot = (double *)calloc(NBin*NBin,sizeof(double));
00122     for(int p=0;p<pNPart();p++){
00123       int v = (int)(pPos(p,CLat1)/pEdge(CLat1)*NBin);
00124       if( v < 0 || v > NBin){ printf("Oi 0 < %d < %d\n",v,NBin);return 1;} 
00125       int vv = (int)(pPos(p,CLat2)/pEdge(CLat2)*NBin);
00126       if( vv < 0 || vv > NBin){ printf("Oi 0 < %d < %d\n",vv,NBin);return 1;} 
00127       Plot[v*NBin+vv] += 1.;
00128     }
00129     for(int v=0;v<NBin;v++)
00130       for(int vv=0,n=0;vv<NBin;vv++)
00131    PlotAv[v*NBin+vv] += Plot[v*NBin+vv];
00132     for(int v=1;v<NBin-1;v++){
00133       for(int vv=1,n=0;vv<NBin-1;vv++){
00134    if(Plot[v*NBin+vv] > Threshold){
00135      NPoint++;
00136      //printf("(%d %d) %d %d %d %d\n",v,vv,Plot[v-1][vv],Plot[v+1][vv],Plot[v][vv-1],Plot[v][vv+1]);
00137      if(Plot[(v-1)*NBin+vv] > Threshold && v > 0)
00138        n++;
00139      if(Plot[(v+1)*NBin+vv] > Threshold && v < NBin -1)
00140        n++;
00141      if(Plot[v*NBin+vv-1] > Threshold && vv > 0)
00142        n++;
00143      if(Plot[v*NBin+vv+1] > Threshold && vv < NBin -1)
00144        n++;
00145      if(n == 4){
00146        Area += 1.;
00147      }
00148      else if(n < 4  && n != 0){
00149        if(n == 3){
00150          Surf += 1.;
00151          Area += 1.;
00152        }
00153        else if(n == 2){
00154          Surf += 2;
00155        }
00156        else if(n == 1){
00157          Surf += sqrt(2);
00158        }    
00159      }
00160      n = 0;
00161    }
00162       }
00163     }
00164     AreaAv += Area;
00165     SurfAv += Surf;
00166     free(Plot);
00167   }
00168   printf("\n");
00169   //-----------normalize-------------------
00170   for(int v=0;v<NBin;v++)
00171     for(int vv=0,n=0;vv<NBin;vv++)
00172       if(PlotAv[v*NBin+vv] > Max)
00173    Max = PlotAv[v*NBin+vv];
00174   //------------write-output-----------------
00175   FILE *FileToWrite = NULL;
00176   FileToWrite = fopen("Surface.xyz","w");
00177   double FNorma = 1./(double)(NFile[1] - NFile[0]);
00178   double LengthConv = pEdge(0)*pEdge(1)/(double)(NBin*NBin);
00179   fprintf(FileToWrite,"#AreaTot %lf Area %lf Surf %lf Threshold %d NChain %d NBin %d\n",pEdge(CLat1)*pEdge(CLat2),AreaAv*FNorma*LengthConv,SurfAv*FNorma*LengthConv,Threshold,pNChain(),NBin);
00180   for(int v=0;v<NBin;v++)
00181     for(int vv=0,n=0;vv<NBin;vv++)
00182       if(PlotAv[v*NBin+vv] > 0.)
00183    fprintf(FileToWrite,"%lf %lf %lf\n",(double)v/NBin*pEdge(CLat1),(double)vv/NBin*pEdge(CLat2),PlotAv[v*NBin+vv]/Max);
00184   fclose(FileToWrite);
00185   free(PlotAv);
00186   return 0;
00187 }
00188 int ElPoly::From3To2d(int NSample,double Param){
00189   char FName[120];
00190   for(int f=NFile[0];f<NFile[1];f++){
00191     Processing(f);
00192     if(OpenRisk(cFile[f],BF_NANO));
00193     sprintf(FName,"ProjOnSurf%05d.dat",f);
00194     FILE *FLine = fopen(FName,"w");
00195     fprintf(FLine,"# l(%lf %lf %lf) d[part]\n",1.,pEdge(CLat2)*pInvEdge(CLat1),.5);
00196     for(int p=0;p<pNPart();p++){
00197       if(pType(p) != 0) continue;
00198       fprintf(FLine,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",pPos(p,CLat1)*pInvEdge(CLat1),pPos(p,CLat2)*pInvEdge(CLat1),.5,pVel(p,0),pVel(p,1),pVel(p,2));
00199     }
00200     fclose(FLine);
00201   }
00202   printf("\n");
00203 }
00204 int ElPoly::From2To1d(int Coord){
00205   if(Coord != 0 && Coord != 1){
00206     printf("Coordinate accepted are 0 or 1 not %d\n",Coord);
00207     return 1;
00208   }
00209   int SecCoord = (Coord+1)%2;
00210   double *Line = (double *) calloc(3*NEdge,sizeof(double));
00211   double *Count = (double *) calloc(3*NEdge,sizeof(double));
00212   for(int p=0;p<pNPart();p++){
00213     int vz = (int)((pPos(p,Coord)+0.001)*pInvEdge(Coord)*NEdge);
00214     if(vz < 0 || vz >= NEdge) continue;
00215     for(int d=0;d<3;d++){
00216       Line[vz*3+d] += pVel(p,d);
00217       Count[vz*3+d] += 1.;
00218     }
00219   }
00220   for(int vz=0;vz<NEdge;vz++){
00221     for(int d=0;d<3;d++){
00222       double Norm = Count[vz*3+d] > 0. ? Count[vz*3+d] : 1.;
00223       //Line[vz*3+d] /= Norm;
00224       Line[vz*3+d] *= pEdge(SecCoord)/(double)NEdge;
00225     }
00226   }
00227   FILE *FLine = fopen("ProjOnLine.dat","w");
00228   for(int v=0;v<NEdge;v++){
00229     fprintf(FLine,"%lf %lf %lf %lf\n",v*pEdge(Coord)/(double)NEdge,Line[v*3+0],Line[v*3+1],Line[v*3+2]);
00230   }
00231   fclose(FLine);
00232   free(Line);
00233   free(Count);
00234 }
00235 int ElPoly::From3To1d(int Coord){
00236 
00237 
00238 }
00239 int ElPoly::RadialShell(int NBin){
00240   double Volume=0;//Global constant
00241   double Area=0.;
00242   double Norm=0.;
00243   double **Plot;
00244   double Ypsilon = 0.;
00245   double InvNFile  = 1./(double)(NFile[1]-NFile[0]);
00246   Plot = (double **)calloc(NBin,sizeof(double));
00247   for(int i=0;i<NBin;i++){
00248     *(Plot+i) = (double *)calloc(NBin,sizeof(double));
00249   }
00250   for(int f=NFile[0];f<NFile[1];f++){
00251     Processing(f);
00252     if(OpenRisk(cFile[f],BF_PART))return 0;
00253     SetEdge(MIN((pEdge(CLat1)*.5),(pEdge(CLat2)*.5)),3);
00254     for(int p=0;p<pNPart();p++){
00255       double Rad = sqrt(SQR((pPos(p,CLat1)-pCm(CLat1))) 
00256          + SQR((pPos(p,CLat2)-pCm(CLat2))) );
00257       int vr = (int)(NBin*Rad*pInvEdge(3));
00258       int vz = (int)(NBin*pPos(p,CNorm)*pInvEdge(CNorm));
00259       //printf("%lf %lf %d %d \n",Rad,pPos(p,CNorm),v,vv);
00260       if( vr < 0 || vr >= NBin) continue;
00261       if( vz < 0 || vz >= NBin) continue;
00262       Plot[vr][vz] += InvNFile;
00263     }
00264     Ypsilon += pEdge(2)-pCm(2)-1.;
00265   }
00266   printf("\n");
00267   FILE *FileToWrite = fopen("RadialShell.xye","w");
00268   //fprintf(FileToWrite,"# l(%.2f %.2f %.2f) v[60] d[chain]\n",Gen->Edge[0],Gen->Edge[1],Gen->Edge[2]);
00269   fprintf(FileToWrite,"%lf %lf %lf\n",0.,0.,0.);
00270   fprintf(FileToWrite,"%lf %lf %lf\n",pEdge(0),pEdge(1),1.);
00271   double Max=0.;
00272   for(int i=0;i<NBin;i++){
00273     for(int j=0;j<NBin;j++){
00274       if(Plot[i][j] > Max) Max = Plot[i][j];
00275     }
00276   }
00277   int th=0;
00278   for(int vz=0;vz<NBin;vz++){
00279     //for(int vr=NBin-1;vr>=0;vr--){
00280     for(int vr=0;vr<NBin;vr++){
00281       if(Plot[vr][vz] > 0.){
00282    fprintf(FileToWrite,"%lf %lf %lf\n",(double)vr/NBin*pEdge(3),(double)vz/NBin*pEdge(2),Plot[vr][vz]/Max);
00283    Norm += Plot[vr][vz];
00284    th++;
00285    if(th > 4){
00286      th = 0;
00287      break;
00288    }
00289       }
00290     }
00291   }
00292   fclose(FileToWrite);
00293   return 0;
00294   Mat->Ypsilon = Ypsilon*InvNFile;
00295   double Vol = 1.;//pNPart()/(Gen->rho/(double)pNPCh()*CUB(5.));
00296   Mat->PreFact = 3./8.*pow((3.*Vol)/DUE_PI,1./3.);
00297   Mat->Func = &Matematica::ContactAngle;
00298   int NRadici = 4;
00299   printf("Volume  %lf Cm %lf Area %lf # to Invert %lf\n",(double)pNPart()/10.,pCm(2),Area,pow(DUE_PI/(2.*pNPart()/10.),.25));
00300   double *Radici;
00301   Radici = (double *)malloc(NRadici*sizeof(double));
00302   int NZeri = Mat->Zeri(0.,DUE_PI/2.,Radici,NRadici);
00303   for(int i=0;i<NZeri;i++){
00304     Radici[i] *= 360./DUE_PI;
00305     printf("Angle %lf\n",Radici[i]);
00306   }
00307   if(NZeri == 0){
00308     printf("The function has no real roots\n");
00309   }
00310   free(Plot);
00311   return 0;
00312 }
00313 void ElPoly::IsoSurf(int NSample,double *IsoLevel,int NIso){
00314   int NPairF = NFile[1]-NFile[0];
00315   double OldPos[3] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2)};
00316   double DensEl = CUB(NSample)/(pVol()*NPairF);
00317   double Dens  = 13.3;
00318   FILE *FNano = fopen("NanoPos.txt","w");
00319   //IsoLevel *= NPairF;
00320   for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
00321     double *Plot = (double *)calloc(CUBE(NSample),sizeof(double));
00322     double Min = 0.;
00323     double Max = 0.;
00324     VAR_TRIANGLE *Triang = NULL;
00325     double Pos[3];
00326     for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
00327       Processing(f);
00328       if(OpenRisk(cFile[f],BF_PART))return;
00329       for(int b=0;b<pNBlock();b++){
00330    for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
00331      if(pType(p) != 0)continue;
00332      for(int d=0;d<3;d++){
00333        Pos[d] = pPos(p,d) - (pNanoPos(0,d)-.5*pEdge(d));
00334        Pos[d] -= floor(Pos[d]*pInvEdge(d))*pEdge(d);
00335      }
00336      int sx = (int)(Pos[0]*pInvEdge(0)*NSample);
00337      int sy = (int)(Pos[1]*pInvEdge(1)*NSample);
00338      int sz = (int)(Pos[2]*pInvEdge(2)*NSample);
00339      int sTot = (sx*NSample+sy)*NSample+sz;
00340      Plot[sTot] += DensEl;
00341      if(Max < Plot[sTot]) Max = Plot[sTot];
00342      if(Min > Plot[sTot]) Min = Plot[sTot];
00343    }
00344       }
00345     }
00346     Matrice Mask(3,3,3);
00347     Mask.FillGaussian(.5,3.);
00348     Mask.Print();
00349     int NDim = 3;
00350     int IfMinImConv = 1;
00351     Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
00352     Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
00353     // ConvoluteMatrix(Plot,NSample,&Mask,3);
00354     // ConvoluteMatrix(Plot,NSample,&Mask,3);
00355     char FString[256];
00356     sprintf(FString,"IsoSurf%05d.dat",ff);
00357     FILE *F2Write = fopen(FString,"w");
00358     fprintf(F2Write,"#l(%lf %lf %lf) v[%d] d[polygon]\n",pEdge(0),pEdge(1),pEdge(2),NSample);
00359     HeaderNano(F2Write);
00360     int NTri = 0;
00361     for(int i=0;i<NIso;i++){
00362       printf("Min %lf Max %lf IsoLevel %lf DensEl %lf\n",Min,Max,IsoLevel[i],DensEl);
00363       Triang = MarchingCubes(Plot,NSample,IsoLevel[i],&NTri);
00364       for(int t=0;t<NTri;t++){
00365    for(int i=0;i<3;i++){
00366      int l1 = t*3 + (i+1)%3;
00367      int l2 = t*3 + (i+2)%3;
00368      for(int d=0;d<3;d++) Pos[d] = Triang[t].p[i].x[d];
00369      int sx = (int)(Pos[0]*pInvEdge(0)*NSample);
00370      int sy = (int)(Pos[1]*pInvEdge(1)*NSample);
00371      int sz = (int)(Pos[2]*pInvEdge(2)*NSample);
00372      int sTot = (sx*NSample+sy)*NSample+sz;
00373      fprintf(F2Write,"{t[%d %d %d] ",sTot,t,0);
00374      fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
00375      fprintf(F2Write,"v(%lf %lf %lf) ",Triang[t].n[i].x[0],Triang[t].n[i].x[1],Triang[t].n[i].x[2]);
00376      fprintf(F2Write,"l[%d] l[%d]}\n",l1,l2);
00377    }
00378       }
00379     }
00380     fclose(F2Write);
00381     free(Plot);
00382     free(Triang);continue;
00383     int NVertex = CUBE(2*NSample-1);
00384     double *WeightL = (double *) calloc(NVertex,sizeof(double));
00385     NormalWeight(Triang,WeightL,NSample,NTri);
00386     double CmStalk[3] = {0.,0.,0.};//OldPos[0],OldPos[1],OldPos[2]};
00387     double CountStalk = 0.;
00388     for(int t=0;t<NTri;t++){
00389       for(int v=0;v<3;v++){
00390    int vRef = Triang[t].v[v];
00391    for(int d=0;d<3;d++){
00392      CmStalk[d] = Triang[t].p[v].x[d]*WeightL[vRef];
00393    }
00394    CountStalk += WeightL[vRef];
00395       }
00396     }
00397     free(WeightL);
00398     free(Triang);
00399     if(CountStalk <= 0.) CountStalk = 1.;
00400     for(int d=0;d<3;d++){
00401       CmStalk[d] /= CountStalk;
00402     }
00403     pPos(CmStalk);
00404     SetNNano(1);
00405     for(int d=0;d<3;d++){
00406       Nano->Pos[d]   = CmStalk[d];
00407       OldPos[d] = Nano->Pos[d];
00408     }
00409     SetNanoBkf(0);
00410     Nano->Axis[0] = 0.;
00411     Nano->Axis[1] = 0.;
00412     Nano->Axis[2] = 1.;
00413     Nano->Rad     = .5;
00414     Nano->Height  = 4.;
00415     Nano->Hamaker = 1.;
00416     Nano->OffSet  = 1.;
00417     Nano->Shape = SHAPE_STALK;    
00418     for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
00419       char String[120];
00420       StringNano(String,0);
00421       fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
00422       fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
00423       //command("chmod u+x %s\n","NanoPos.txt");
00424       //SubNanoHeader(cFile[f]);
00425     }
00426     printf("\n");
00427   }
00428   fclose(FNano);
00429 }
00430 void ElPoly::IsoLine(int NSample,double *IsoLevel,int NIso,int How){
00431   int NPairF = 1;//NFile[1]-NFile[0];
00432   double OldPos[3] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2)};
00433   double DensEl = SQR(NSample)/(pVol()*NPairF);
00434   double Dens  = 13.3;
00435   //IsoLevel *= NPairF;
00436   for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
00437     double *Plot = (double *)calloc(SQR(NSample),sizeof(double));
00438     double Min = 0.;
00439     double Max = 0.;
00440     double Pos[3];
00441     for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
00442       Processing(f);
00443       if(OpenRisk(cFile[f],BF_PART))return;
00444       for(int b=0;b<pNBlock();b++){
00445    for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
00446      if(pType(p) != 0)continue;
00447      for(int d=0;d<3;d++){
00448        Pos[d] = pPos(p,d) - (pNanoPos(0,d)-.5*pEdge(d));
00449        Pos[d] -= floor(Pos[d]*pInvEdge(d))*pEdge(d);
00450      }
00451      int sx = (int)(Pos[CLat1]*pInvEdge(CLat1)*NSample);
00452      int sy = (int)(Pos[CLat2]*pInvEdge(CLat2)*NSample);
00453      int sTot = sx*NSample+sy;
00454      if(How==0)//particle file
00455        Plot[sTot] += DensEl;
00456      else
00457        Plot[sTot] += pPos(p,CNorm);
00458      if(Max < Plot[sTot]) Max = Plot[sTot];
00459      if(Min > Plot[sTot]) Min = Plot[sTot];
00460    }
00461       }
00462       printf("Min %lf Max %lf DensEl %lf\n",Min,Max,DensEl);
00463     }
00464     Matrice Mask(5,5);
00465     Mask.FillGaussian(.5,3.);
00466     Mask.Print();
00467     int NDim = 2;
00468     int IfMinImConv = 1;
00469     Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
00470     Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
00471     char FString[60];
00472     sprintf(FString,"IsoSurf%05d.dat",ff);
00473     FILE *F2Write = fopen(FString,"w");
00474     fprintf(F2Write,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(0),pEdge(1),pEdge(2),NSample);
00475     HeaderNano(F2Write);
00476     IsoLine(F2Write,Plot,NSample,IsoLevel,NIso);
00477     free(Plot);
00478     fclose(F2Write);
00479   }
00480 }
00481 void ElPoly::IsoLine(FILE *F2Write,double *Plot,int NSample,double *IsoLevel,int NIso){
00482   VAR_LINE *Triang = NULL;
00483   int NTri = 0;
00484   int p=0;
00485   char FName[60];
00486   double Pos[3];
00487   for(int i=0;i<NIso;i++){
00488     Triang = MarchingSquares(Plot,NSample,IsoLevel[i],&NTri);
00489     ConnectLineChain(Triang,NSample,NTri);
00490     sprintf(FName,"IsoLineChain%d.dat",i);
00491     sprintf(cWhat2Draw,"part");
00492     Write(FName);
00493     // for(int c=0;c<pNChain();c++){
00494     //   sprintf(FString,"Line%05dChain%02d.dat",ff,c);
00495     //   FILE *FLine = fopen(FString,"w");
00496     // while(
00497     // for(int p=0;p<pNChain();p++){
00498     //   if
00499     //   fprintf(FString,"%lf %lf\n"
00500     //   fclose(FLine);
00501     // }
00502     if(1==1){//Isoline without ordering
00503       for(int t=0;t<NTri;t++){
00504    for(int v=0;v<2;v++){
00505      for(int d=0;d<3;d++) Pos[d] = Triang[t].p[v].x[d];
00506      int sTot = Triang[t].v[v];
00507      fprintf(F2Write,"{t[%d %d %d] ",sTot,t,0);
00508      fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
00509      if(v==0)
00510        fprintf(F2Write,"l[%d]}\n",p+1);
00511      else
00512        fprintf(F2Write,"\n");
00513      p++;
00514    }
00515       }
00516     }
00517     if(1==0){//Isoline separated between inside and outside an elips
00518       double Radx = 20.;
00519       double Rady = 12.;
00520       double Center[3] = {0.,.5*pEdge(CLat2),0.};
00521       int NIt = 100;
00522       double DistElips = SQR(Radx)+SQR(Rady);
00523       for(int i=0;i<NIt;i++){
00524    double Ang = i/(double)NIt*M_PI - M_PI*.5;
00525    double x = Center[0] + Radx*cos(Ang);
00526    double y = Center[1] + Rady*sin(Ang);
00527    fprintf(F2Write,"{t[%d %d %d] ",p++,2,2);
00528    fprintf(F2Write,"x(%lf %lf %lf) }\n",x,y,0.);     
00529       }
00530       for(int t=0;t<NTri;t++){
00531    for(int v=0;v<2;v++){
00532      for(int d=0;d<3;d++) Pos[d] = Triang[t].p[v].x[d];
00533      int sTot = Triang[t].v[v];
00534      int type = 0;
00535      if(fabs(Pos[0] - Center[0]) < Radx && fabs(Pos[1] - Center[1]) < Rady) type = 1;
00536      fprintf(F2Write,"{t[%d %d %d] ",sTot,t,type);
00537      fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
00538      if(v==0)
00539        fprintf(F2Write,"l[%d]}\n",p+1);
00540      else
00541        fprintf(F2Write,"\n");
00542      p++;
00543    }
00544       }
00545     }
00546   }
00547   if(1==1){//to write the continuum density values
00548     for(int sx=0;sx<NSample;sx++){
00549       for(int sy=0;sy<NSample;sy++){
00550    double x = pEdge(CLat1)*sx/(double)NSample;
00551    double y = pEdge(CLat2)*sy/(double)NSample;
00552    int sTot = sx*NSample+sy;
00553    fprintf(F2Write,"{t[%d %d %d] ",sTot,NTri+sTot/2,1);
00554    fprintf(F2Write,"x(%lf %lf %lf) }\n",x,y,Plot[sx*NSample+sy]);
00555       }
00556     }
00557   }
00558   if(1==0){//Write the chains separately
00559     for(int c=0;c<pNChain();c++){
00560       sprintf(FName,"Chain%d.dat",c);
00561       FILE *FChain = fopen(FName,"w");
00562       for(int p=0;p<pNPart();p++){
00563    if(Pm[p].CId != c) continue;
00564    fprintf(FChain,"%lf %lf\n",Pm[p].Pos[0],Pm[p].Pos[1]);
00565       }
00566     }
00567   }
00568   free(Triang);
00569 }
00570 void ElPoly::FetchPore(){
00571   FILE *FNano = fopen("NanoPos.sh","w");
00572   double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
00573   FILE *StalkArea = fopen("PoreArea.dat","w");
00574   fprintf(StalkArea,"#time rad asy\n");
00575   for(int f=NFile[0];f<NFile[1];f++){
00576     Processing(f);
00577     if(OpenRisk(cFile[f],BF_PART)) return;
00578     double Asy = PorePos();
00579     Nano[0].Shape = ShapeId("pore");
00580     Nano->Height = .2;
00581     char String[120];
00582     StringNano(String,0);
00583     fprintf(StalkArea,"%lf %lf %lf\n",pTime(),Nano->Rad,Asy);
00584     fprintf(FNano,"sed '/pore/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
00585     fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
00586     //command("chmod u+x %s\n","NanoPos.txt");
00587     //SubNanoHeader(cFile[f]);
00588   }
00589   fclose(FNano);
00590   fclose(StalkArea);
00591   printf("\n");
00592 }
00596 void ElPoly::FetchStalk(){
00597   FILE *FNano = fopen("NanoPos.sh","w");
00598   double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
00599   FILE *StalkArea = fopen("StalkArea.dat","w");
00600   char FName[60];
00601   for(int f=NFile[0];f<NFile[1];f++){
00602     Processing(f);
00603     if(OpenRisk(cFile[f],BF_NO)) return;
00604     SetNNano(1);
00605     if(StalkPos(OldPos)) continue;
00606     char String[120];
00607     StringNano(String,0);
00608     fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
00609     fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
00610     fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
00611     sprintf(FName,"Centered%05d.dat",f);
00612     //Write(FName);
00613     //command("chmod u+x %s\n","NanoPos.txt");
00614     //SubNanoHeader(cFile[f]);
00615   }
00616   fclose(FNano);
00617   fclose(StalkArea);
00618   printf("\n");
00619 }
00624 void ElPoly::StalkArea(){
00625   FILE *FNano = fopen("NanoPosA.sh","w");
00626   FILE *StalkArea = fopen("StalkArea.dat","w");
00627   FILE *AreaS = fopen("AreaS.dat","w");
00628   double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
00629   int NBin = 36;
00630   double *Count = (double *)calloc(NBin*NBin,sizeof(double));
00631   double RadTorus = 1.;//Nano->Rad;
00632   //fprintf(AreaS,"#l(%lf %lf %lf) \n",2.*Nano->Height,2.*Nano->Height,10.);
00633   fprintf(AreaS,"#l(%lf %lf %lf) \n",pEdge(0),pEdge(1),pEdge(2));
00634   HeaderNano(AreaS);
00635   char FName[60];
00636   for(int f=NFile[0];f<NFile[1];f++){
00637     Processing(f);
00638     if(OpenRisk(cFile[f],BF_NO)) return;
00639     SetNNano(1);
00640     if(StalkPos(OldPos)) continue;
00641     //Nano->Pos[CNorm] = pCm(CNorm);
00642     //Nano->Pos[CNorm] -= floor(Nano->Pos[CNorm]*pInvEdge(CNorm))*pEdge(CNorm);
00643     double Cm[3] = {0.,0.,0.};
00644     double CountCm = 0;
00645     int nPart = 0;
00646     double Pos[3];
00647     //counts the particles inside the torus
00648     for(int p=0;p<pNPart();p++){
00649       for(int d=0;d<3;d++){
00650    Pos[d] = pPos(p,d) - Nano->Pos[d];
00651    Pos[d] -= floor(pPos(p,d)*pInvEdge(d))*pEdge(d);
00652       }
00653       double Rad = hypot(Pos[CLat1],Pos[CLat2]);
00654       if(Rad > Nano->Height) continue;
00655       if(fabs(Pos[CNorm]) > RadTorus) continue;
00656       // fprintf(AreaS,"{t[%d 0 %d] x(%lf %lf %lf)",nPart++,pType(p),pPos(p,0),pPos(p,1),pPos(p,2));
00657       // if(Ln[p].NLink > 0) fprintf(AreaS,"l[%d]",p-Ln[p].Link[0]+nPart+1);
00658       // fprintf(AreaS,"}\n");
00659       if(pType(p) == 1) continue;
00660       for(int d=0;d<3;d++){
00661    Cm[d] += pPos(p,d);
00662       }
00663       CountCm += 1.;
00664       int vx = (int)(Pos[CLat1]/(2.*Nano->Height)*NBin);
00665       vx += NBin/2;
00666       if(vx < 0 || vx >= NBin) continue;
00667       int vy = (int)(Pos[CLat2]/(2.*Nano->Height)*NBin);
00668       vy += NBin/2;
00669       if(vy < 0 || vy >= NBin) continue;
00670       Count[vx*NBin+vy] += 1.;
00671     }
00672     double Area = 0.;
00673     double Norm = 0.;
00674     for(int vx=0;vx<NBin;vx++){
00675       for(int vy=0;vy<NBin;vy++){
00676       if(Count[vx*NBin+vy] < 1.) continue;
00677    // double x = vx*Nano->Height*2./(double)NBin;
00678    // double y = vy*Nano->Height*2./(double)NBin;
00679    // fprintf(AreaS,"{t[%d 0 2] x(%lf %lf %lf)}\n",nPart++,x+pNanoPos(0,0)-Nano->Height,y+pNanoPos(0,1)-Nano->Height,pNanoPos(0,2));
00680       Area += 1.;
00681       }
00682     }
00683     if(CountCm <= 0.){
00684       printf("No particles in the torus %s\n",cFile[f]);
00685       return;
00686     }
00687     Nano->Area = SQR(2.*Nano->Height)*Area/(double)(SQR(NBin));
00688     for(int d=0;d<3;d++){
00689       Cm[d] /= CountCm;
00690     }
00691     Cm[CNorm] = pCm(CNorm);
00692     //fprintf(AreaS,"%lf %lf %lf\n",Cm[0]-Nano->Pos[0],Cm[1]-Nano->Pos[1],Cm[2]-Nano->Pos[2]);
00693     for(int d=0;d<3;d++){
00694       Nano->Pos[d] = Cm[d] - floor(Cm[d]*pInvEdge(d))*pEdge(d);
00695       fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
00696     }
00697     SetNanoBkf(0);
00698     Nano->Height = Nano->Rad + sqrt(Nano->Area/DUE_PI);
00699     char String[120];
00700     StringNano(String,0);
00701     fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
00702     fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
00703     fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
00704     sprintf(FName,"Centered%05d.dat",f);
00705     //Write(FName);
00706     //HeaderNano(AreaS);
00707   }
00708   fclose(AreaS);
00709   fclose(StalkArea);
00710   fclose(FNano);
00711   free(Count);
00712   printf("\n");
00713 }
00714 void ElPoly::AvSnap(){
00715   int NPairF = 10;
00716   double InvAv = 1./(double)NPairF;
00717   double Norm = 1./(double)(NFile[1]-NFile[0]);
00718   char FName[256];
00719   // for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
00720   PART *Pn = (PART *)calloc(pNPart(),sizeof(PART));
00721   //for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
00722   for(int f=NFile[0];f<NFile[1];f++){
00723     Processing(f);
00724     if(OpenRisk(cFile[f],BF_NO))return;
00725     for(int p=0;p<pNPart();p++){
00726       for(int d=0;d<3;d++){
00727    Pn[p].Pos[d] += pPos(p,d);
00728       }
00729     }
00730   }
00731   for(int p=0;p<pNPart();p++){
00732     double Pos[3] = {Pn[p].Pos[0]*Norm,Pn[p].Pos[1]*Norm,Pn[p].Pos[2]*Norm};
00733     SetPos(p,Pos);
00734   }
00735   printf("\n");
00736   free(Pn);
00737   sprintf(FName,"Av%s",cFile[NFile[0]]);
00738   Write(FName);
00739 }
00740 void ElPoly::Sample(int NSample){
00741   MOMENTI m1 = SampleSurfaceMem(NSample);
00742   FILE *FWrite = fopen("SolSim.dat","w");
00743   double InvNSample = 1./(double)NSample;
00744   fprintf(FWrite,"# l(%.2f %.2f %.2f) d[part]\n",pEdge(0),pEdge(1),pEdge(2));
00745   for(int sx=0;sx<NSample;sx++){
00746     int sy = 0;
00747     PlotMem[sx*NSample + sy] = m1.Uno;
00748     PlotMem[sy*NSample + sx] = m1.Uno;
00749     sy = NSample - 1;
00750     PlotMem[sx*NSample + sy] = m1.Uno;
00751     PlotMem[sy*NSample + sx] = m1.Uno;
00752   }
00753   printf("sample average: %lf\n",m1.Uno);
00754   double ConvFact = 1.45/m1.Uno;
00755   for(int sx=0;sx<NSample;sx++){
00756     double x = sx*InvNSample*pEdge(CLat1);
00757     int sx1 = sx + 1 == NSample ? 0 : sx + 1;
00758     for(int sy=0;sy<NSample;sy++){
00759       double y = sy*InvNSample*pEdge(CLat2);
00760       int sy1 = sy + 1 == NSample ? 0 : sy + 1;
00761       int p = sx*NSample+sy;
00762       int l1 = sx1*NSample+sy + SQR(NSample);
00763       int l2 = sx*NSample+sy1 + SQR(NSample);
00764       //fprintf(FWrite,"{t[%d 0 2] x(%lf %lf %lf) l[%d] l[%d]}\n",p,x,y,PlotMem[p],sx1*NSample+sy,sx*NSample+sy1);
00765       fprintf(FWrite,"{t[%d 0 3] x(%lf %lf %lf) l[%d] l[%d]}\n",p,x,y,(PlotMem[p]-m1.Uno)*ConvFact,l1,l2);
00766       //fprintf(FWrite,"{t[%d 0 2] x(%lf %lf %lf)}\n",p,x,y,PlotMem[p]-m1.Uno);
00767     }
00768   }
00769 }