Allink: src/ElPoly/ElPolyOutput.cpp Source File

00001 #include "ElPoly.h"
00002 
00003 void ElPoly::Conv2Tecplot(int NBin,int How){
00004   int NType = 3;
00005   double Round = 0.001;
00006   int Nx = MIN(NBin,NEdge);
00007   int Ny = MIN(NBin,NEdge);
00008   double **Plot = (double **)calloc(NType,sizeof(double));
00009   double **Count = (double **)calloc(NType,sizeof(double));
00010   for(int t=0;t<NType;t++){
00011     Plot[t] = (double *)calloc(Nx*Ny,sizeof(double));
00012     Count[t] = (double *)calloc(Nx*Ny,sizeof(double));    
00013   }
00014   FILE *TecPlot = fopen("TecPlot.dat","w");
00015   if(How == 0){//Dens
00016     fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"oil density\",\"density phob\",\"density phil\"\n");
00017     fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
00018   }
00019   else if(How == 1){
00020     fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"pressure\",\"density phob\",\"density phil\"\n");
00021     fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
00022   }
00023   else if(How == 2){
00024     fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"height\",\"density phob\",\"density phil\"\n");
00025     fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
00026   }
00027   for(int p=0;p<pNPart();p++){
00028     int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*Nx);
00029     if(vx < 0 || vx >= Nx) continue;
00030     int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*Ny);
00031     if(vy < 0 || vy >= Ny) continue;
00032     if(How == 0){//Dens
00033       int t = Pm[p].Typ;
00034       t = (t+1)%3;
00035       // if(t==2){t=0;Plot[1][vx*Nx+vy] += -1.;}
00036       // else 
00037    Plot[t][vx*Nx+vy] += pPos(p,2);
00038       Count[t][vx*Nx+vy] += 1.;
00039     }
00040     else if(How == 1){
00041       for(int t=0;t<3;t++){
00042    Plot[t][vx*Nx+vy] += pVel(p,t);
00043    Count[t][vx*Nx+vy] += 1.;
00044       }
00045     }
00046     else if(How == 2){
00047       int t = 0;//Pm[p].Typ;
00048       Plot[t][vx*Nx+vy] += pPos(p,2);
00049       Count[t][vx*Nx+vy] += 1.;
00050     }
00051   }
00052   Matrice Mask(5,5);
00053   Mask.FillGaussian(.5,3.);
00054   int NDim = 2;
00055   int IfMinImConv = 0;
00056   for(int t=0;t<3;t++){
00057     //Mask.ConvoluteMatrix(Plot[t],Nx,2);
00058     Mask.ConvoluteMatrix(Plot[t],Nx,NDim,IfMinImConv);
00059   }
00060   for(int vx=0;vx<Nx;vx++){
00061     for(int vy=0;vy<Ny;vy++){
00062       double r = vx*pEdge(0)/(double)Nx;
00063       double z = vy*pEdge(1)/(double)Ny - pEdge(1)*.5;
00064       double Norm0 = Count[0][vx*Nx+vy] > 0. ? 1./Count[0][vx*Nx+vy] : 1.;
00065       double Norm1 = Count[1][vx*Nx+vy] > 0. ? 1./Count[1][vx*Nx+vy] : 1.;
00066       double Norm2 = Count[2][vx*Nx+vy] > 0. ? 1./Count[2][vx*Nx+vy] : 1.;
00067       fprintf(TecPlot,"%lf %lf %lf %lf %lf\n",r,z,
00068          Plot[0][vx*Nx+vy]*Norm0,Plot[1][vx*Nx+vy]*Norm1,Plot[2][vx*Nx+vy]*Norm2);
00069     }
00070   }
00071   // if(Pm[p].Vel[0]+Pm[p].Vel[1]+Pm[p].Vel[2]<.1){
00072   //   Dens1 = -5.;
00073   //   Dens2 = -5.;
00074   //   Dens3 = -5.;
00075   // }
00076   // if(pPos(p,0) > 0 && pPos(p,0) < 5.)
00077   //   if(pPos(p,1)>13. && pPos(p,1)<19.)
00078   //  if(Pm[p].Vel[0] < 0.01){
00079   //    Dens1 = -3.;
00080   //    Dens2 = -3.;
00081   //    Dens3 = -3.;
00082   //  }
00083   fclose(TecPlot);
00084   for(int t=0;t<NType;t++){
00085     free(Plot[t]);
00086     free(Count[t]);
00087   }
00088   free(Plot);
00089   free(Count);
00090 }
00091 void ElPoly::Conv2Vmd(){
00092   FILE *OutVmd = fopen("Sim.vtf","w");
00093   int cOff = 0;
00094   int pOff = 0;
00095   char Shape[30];
00096   int NanoPoint = 10;
00097   for(int b=0;b<pNBlock();b++){
00098     int bType = 0;
00099     if(!strncmp(Block[b].Name,"PEP",3)){
00100       bType = 2;
00101     }    
00102     for(int c=cOff;c<cOff+pNChain(b);c++){
00103       for(int p=pOff;p<pOff+pNPCh(b);p++){
00104    int Type = pType(p) + bType;
00105    if(Type == 0)
00106      fprintf(OutVmd,"a %d r 0.8 n A resid %d res %s\n",p,pChain(p),Block[b].Name);
00107    else if(Type == 1)
00108      fprintf(OutVmd,"a %d r 0.8 n B resid %d res %s\n",p,pChain(p),Block[b].Name);
00109    else if(Type == 2)
00110      fprintf(OutVmd,"a %d r 0.8 n D resid %d res %s\n",p,pChain(p),Block[b].Name);
00111    else if(Type == 3)
00112      fprintf(OutVmd,"a %d r 0.8 n E resid %d res %s\n",p,pChain(p),Block[b].Name);
00113       }
00114       pOff += pNPCh(b);
00115       fprintf(OutVmd,"b %d::%d\n",c*pNPCh(b),(c+1)*pNPCh(b)-1);
00116     }
00117     cOff += pNChain(b);
00118   }
00119   for(int n=0;n<pNNano();n++){
00120     ShapeId(Nano[n].Shape,Shape);
00121     for(int i=0;i<NanoPoint;i++){
00122       fprintf(OutVmd,"a %d r %lf n Nano resid %d res %s\n",n*NanoPoint+pNPart()+i,Nano[n].Rad,cOff+n,Shape);
00123     }
00124     fprintf(OutVmd,"b %d::%d\n",n*NanoPoint+pNPart(),n*NanoPoint+pNPart()+NanoPoint-1);
00125   }
00126   for(int f=NFile[0];f<NFile[1];f++){
00127     Processing(f);
00128     if(OpenRisk(cFile[f],BF_CHAIN)) return ;
00129     fprintf(OutVmd,"\ntimestep\npbc %lf %lf %lf\n",pEdge(0),pEdge(1),pEdge(2));
00130     for(int p=0;p<pNPart();p++){
00131       fprintf(OutVmd,"%lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2));
00132     }
00133     for(int n=0;n<pNNano();n++){
00134       double Pos[3];
00135       for(int d=0;d<3;d++){
00136    Pos[d] = pNanoPos(n,d) - .5*Nano[n].Height*Nano[n].Axis[d];
00137       }
00138       for(int i=0;i<NanoPoint;i++){
00139    for(int d=0;d<3;d++){
00140      Pos[d] += Nano[n].Height/(double)NanoPoint*Nano[n].Axis[d];
00141    }
00142    fprintf(OutVmd,"%lf %lf %lf\n",Pos[0],Pos[1],Pos[2]);
00143       }
00144     }
00145   }
00146 }
00147 //-----------------------POVRAY--------------------------
00148 void ElPoly::DrPartPovRay(int p){
00149   double Pos[3];
00150   for(int d=0;d<3;d++) Pos[d] = pPos(p,d)*InvScaleUn;
00151   int Typ = pType(p) < 6 ? pType(p) : 5;
00152   int Chc = Ch[Pm[p].CId].Type;
00153   if(pType(p) == 0 && CHAIN_IF_TYPE(Chc,CHAIN_ADDED) )Typ = 3;
00154   fprintf(DrawOutFile,"ElSphere(<%.4f, %.4f, %.4f>,",Pos[0],Pos[1],Pos[2]);
00155   fprintf(DrawOutFile,"SphRad,rgbt<%.4f, %.4f, %.4f, %.4f>)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
00156 }
00157 void ElPoly::DrNanoPovRay(int n){
00158   if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_CYL)){
00159     double PosP[3];
00160     double PosN[3];
00161     int Typ = 2;
00162     for(int d=0;d<3;d++){
00163       PosP[d] = pNanoPos(n,d) + .5*Nano[n].Height*Nano[n].Axis[d];
00164       PosN[d] = pNanoPos(n,d) - .5*Nano[n].Height*Nano[n].Axis[d];
00165       PosP[d] *= InvScaleUn;
00166       PosN[d] *= InvScaleUn;
00167     }
00168     fprintf(DrawOutFile,"ElCylinder(<%.2f, %.2f, %.2f>,",PosP[0],PosP[1],PosP[2]);
00169     fprintf(DrawOutFile,"<%.2f, %.2f, %.2f>%.2f,",PosN[0],PosN[1],PosN[2],Nano[n].Rad*InvScaleUn);
00170     fprintf(DrawOutFile,"rgbt<%.4f, %.4f, %.4f, %.4f>,1)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
00171   }
00172   else{
00173     // Point2Shape(Nano[n].Shape);
00174     // DrField(128,SQR(Nano[n].Rad),n,DrawOutFile);
00175   }
00176 }
00177 void ElPoly::DrBondPovRay(double *Pos1,double *Pos2,float *Color){
00178   fprintf(DrawOutFile,"ElCylinder(<%lf, %lf, %lf>,",Pos1[0],Pos1[1],Pos1[2]);
00179   fprintf(DrawOutFile,"<%lf, %lf, %lf>",Pos2[0],Pos2[1],Pos2[2]);
00180   fprintf(DrawOutFile,"CylRad,rgbt<%.4f, %.4f, %.4f, %.4f>,1)\n",Color[0],Color[1],Color[2],Color[3]-1.);
00181 }
00182 void ElPoly::Conv2Povray(){
00183   //particles, lines
00184   HeaderPovRay();
00185   for(int f=NFile[0];f<NFile[1];f++){
00186     SigErr(DrawOutFile != NULL,"DrawOutFile already allocated, can't use the file");
00187     char FName[60];
00188     sprintf(FName,"PovSnap%05d.pov",pStep());
00189     DrawOutFile = fopen(FName,"w");
00190     int ImSize[2] = {1000,1000};
00191     fprintf(DrawOutFile,"// POV 3.x input script : plot.pov\n// command: povray +W%d +H%d -I%s -O%s.tga +P +X +A +FT +C\n",ImSize[0],ImSize[1],FName,FName);
00192     fprintf(DrawOutFile,"#if (version < 3.5)\n#error \"POV3DisplayDevice has been compiled for POV-Ray 3.5 or above.\\nPlease upgrade POV-Ray.\"\n#end\n");
00193     fprintf(DrawOutFile,"#include \"PovHeader.inc\"\n");
00194     fprintf(DrawOutFile,"// System \n");
00195     for(int b=0,NPep=0;b<pNBlock();b++){
00196       for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
00197    DrPartPovRay(p);
00198    int Typ = Pm[p].Typ;
00199    double Pos1[3];
00200    double Pos2[3];
00201    for(int l=0;l<Ln[p].NLink;l++){
00202      int link = Ln[p].Link[l];
00203      if(p == link) continue;
00204      for(int d=0;d<3;d++){
00205        Pos1[d] = pPos(p,d)*InvScaleUn;
00206        Pos2[d] = pPos(p,link)*InvScaleUn;
00207        if(Pos1[d] - Pos2[d] > .5*InvScaleUn)
00208          Pos2[d] += pEdge(d)*InvScaleUn;
00209        else if(Pos1[d] - Pos2[d] < -.5*InvScaleUn)
00210          Pos2[d] -= pEdge(d)*InvScaleUn;
00211      }
00212      DrBondPovRay(Pos1,Pos2,ColorType[Typ]);
00213    }
00214       }
00215     }
00216     for(int n=0;n<pNNano();n++){
00217       DrNanoPovRay(n);
00218     }
00219     fclose(DrawOutFile);
00220   }
00221 }
00222 #include "ElPolyDrawSurf.h"
00223 void ElPoly::DrField(int NGrid,double IsoLevel,int nNano,FILE *FWrite){
00224   int Typ = 2;
00225   double CubeDist[8];
00226   double EdgeVertex[12][3];
00227   double EdgeNormal[12][3];
00228   double InvNGrid = 1./(double)NGrid;
00229   double Pos[3];
00230   double Pos1[3];
00231   double Cm[3];
00232   for(int d=0;d<3;d++){
00233     Cm[d] = .5*pEdge(d)*InvScaleUn;
00234   }
00235   for(int gx=0;gx<NGrid;gx++){
00236     Pos[0] = gx*InvNGrid*pEdge(0);
00237     for(int gy=0;gy<NGrid;gy++){
00238       Pos[1] = gy*InvNGrid*pEdge(1);
00239       for(int gz=0;gz<NGrid;gz++){
00240    Pos[2] = gz*InvNGrid*pEdge(2);
00241    for(int v=0;v<8;v++){
00242      for(int d=0;d<3;d++){
00243        Pos1[d] = Pos[d] + VertCube[v][d]*InvNGrid*pEdge(d);
00244      }
00245      CubeDist[v] = NanoDist2(Pos1,nNano);
00246    }
00247    int Flag = 0;
00248    for(int v=0;v<8;v++){
00249      if(CubeDist[v] <= IsoLevel)
00250        Flag |= 1<<v;
00251    }
00252    int CubeShape = CubeTop[Flag];
00253    if(CubeShape==0) continue;
00254    for(int e=0;e<12;e++){
00255      if(CubeShape & (1<<e)){
00256        double Delta = CubeDist[EdgeConn[e][1]] - CubeDist[EdgeConn[e][0]];
00257        double OffSet = (IsoLevel-CubeDist[EdgeConn[e][0]])/Delta;
00258        if(Delta == 0.0){
00259          OffSet = .5;
00260        }
00261        EdgeNormal[e][0] = NanoDist2(Pos[0]-0.01,Pos[1],Pos[2],nNano) 
00262          - NanoDist2(Pos[0]+0.01,Pos[1],Pos[2],nNano);
00263        EdgeNormal[e][1] = NanoDist2(Pos[0],Pos[1]-0.01,Pos[2],nNano) 
00264          - NanoDist2(Pos[0],Pos[1]+0.01,Pos[2],nNano);
00265        EdgeNormal[e][2] = NanoDist2(Pos[0],Pos[1],Pos[2]-0.01,nNano) 
00266          - NanoDist2(Pos[0],Pos[1],Pos[2]+0.01,nNano);
00267        double Norm = sqrt(SQR(EdgeNormal[e][0])+SQR(EdgeNormal[e][1])+SQR(EdgeNormal[e][2]));
00268        for(int d=0;d<3;d++){
00269          EdgeVertex[e][d] = Pos[d] + (VertCube[EdgeConn[e][0]][d]+OffSet*EdgeDir[e][d])*InvNGrid*pEdge(d);
00270          EdgeVertex[e][d] *= InvScaleUn;
00271          EdgeNormal[e][d] *= 1./Norm;
00272        }
00273      }
00274    }
00275    for(int t=0;t<5;t++){
00276      if(TrConnTable[Flag][3*t] < 0.) break;
00277      fprintf(FWrite,"ElTriangle (");
00278      for(int d=0;d<3;d++){
00279        int v = TrConnTable[Flag][3*t+d];
00280        fprintf(FWrite,"<%.2f, %.2f, %.2f>,",EdgeVertex[v][0]-Cm[0],EdgeVertex[v][1]-Cm[1],EdgeVertex[v][2]-Cm[2]);
00281        fprintf(FWrite,"<%.2f, %.2f, %.2f>,",EdgeNormal[v][0],EdgeNormal[v][1],EdgeNormal[v][2]);
00282      }
00283      fprintf(FWrite,"rgbt<%.3f, %.3f, %.3f, %.3f>)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
00284    }
00285       }
00286     }
00287   }
00288 }
00289 void ElPoly::Conv2rzd(int NBin){
00290   // FILE *FOut = fopen("data.dat","w");
00291   // for(int p=0;p<pNPart();p++){
00292   //   fprintf(FOut,"%lf %lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0));
00293   // }
00294   // fclose(FOut);
00295   // return;
00296   int NType = 3;
00297   double Round = 0.001;
00298   double **Plot = (double **)calloc(NType,sizeof(double));
00299   for(int t=0;t<NType;t++)
00300     Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
00301   double *Count = (double *)calloc(NBin*NBin*NType,sizeof(double));
00302   FILE *TecPlot = fopen("data.dat","w");
00303   for(int p=0;p<pNPart();p++){
00304     int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*NBin);
00305     if(vx < 0 || vx >= NBin) continue;
00306     int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*NBin);
00307     if(vy < 0 || vy >= NBin) continue;
00308     for(int t=0;t<3;t++){
00309       Plot[t][vx*NBin+vy] += pVel(p,t);
00310       Count[(vx*NBin+vy)*NType+t] += 1.;
00311     }
00312   }
00313   //smooth and write dens
00314   Matrice Mask(5,5);
00315   Mask.FillGaussian(.5,3.);
00316   Mask.Print();
00317   int NDim = 2;
00318   int IfMinImConv = 0;
00319   for(int t=0;t<3;t++){
00320     Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
00321     Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
00322   }
00323   for(int vx=0;vx<NBin;vx++){
00324     for(int vy=0;vy<NBin;vy++){
00325       double r = vx*pEdge(0)/(double)NBin;
00326       double z = vy*pEdge(1)/(double)NBin - pEdge(1)*.5;
00327       double Norm0 = Count[(vx*NBin+vy)*NType+0] > 0. ? 1./Count[(vx*NBin+vy)*NType+0] : 1.;
00328       double Norm1 = Count[(vx*NBin+vy)*NType+1] > 0. ? 1./Count[(vx*NBin+vy)*NType+1] : 1.;
00329       double Norm2 = Count[(vx*NBin+vy)*NType+2] > 0. ? 1./Count[(vx*NBin+vy)*NType+2] : 1.;
00330       fprintf(TecPlot,"%lf %lf %lf %lf %lf\n",r,z,
00331          Plot[0][vx*NBin+vy]*Norm0,Plot[1][vx*NBin+vy]*Norm1,Plot[2][vx*NBin+vy]*Norm2);
00332       // fprintf(TecPlot,"%lf %lf %lf \n",r,z,Plot[1][vx*NBin+vy]*Norm0);
00333     }
00334   }
00335   fclose(TecPlot);
00336   for(int t=0;t<NType;t++)
00337     free(Plot[t]);
00338   free(Plot);
00339   free(Count);
00340 }
00341 void ElPoly::Conv2xyzd(int NBin){
00342   // FILE *FOut = fopen("data.dat","w");
00343   // for(int p=0;p<pNPart();p++){
00344   //   fprintf(FOut,"%lf %lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0));
00345   // }
00346   // fclose(FOut);
00347   // return;
00348   int NType = 3;
00349   double Round = 0.001;
00350   double **Plot = (double **)calloc(NType,sizeof(double));
00351   for(int t=0;t<NType;t++)
00352     Plot[t] = (double *)calloc(NBin*NBin*NBin,sizeof(double));
00353   double *Count = (double *)calloc(NBin*NBin*NBin*NType,sizeof(double));
00354   FILE *TecPlot = fopen("data.dat","w");
00355   for(int p=0;p<pNPart();p++){
00356     int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*NBin);
00357     if(vx < 0 || vx >= NBin) continue;
00358     int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*NBin);
00359     if(vy < 0 || vy >= NBin) continue;
00360     int vz = (int)((pPos(p,2)+Round)*pInvEdge(2)*NBin);
00361     if(vz < 0 || vz >= NBin) continue;
00362     for(int t=0;t<3;t++){
00363       Plot[t][(vx*NBin+vy)*NBin+vz] += pVel(p,t);
00364       Count[((vx*NBin+vy)*NBin+vz)*NType+t] += 1.;
00365     }
00366   }
00367   //smooth and write dens
00368   for(int vx=0;vx<NBin;vx++){
00369     for(int vy=0;vy<NBin;vy++){
00370       for(int vz=0;vz<NBin;vz++){
00371    double x = vx*pEdge(0)/(double)NBin;
00372    double y = vy*pEdge(1)/(double)NBin;
00373    double z = vz*pEdge(1)/(double)NBin;
00374    double Norm0 = Count[((vx*NBin+vy)*NBin+vz)*NType+0] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+0] : 1.;
00375    double Norm1 = Count[((vx*NBin+vy)*NBin+vz)*NType+1] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+1] : 1.;
00376    double Norm2 = Count[((vx*NBin+vy)*NBin+vz)*NType+2] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+2] : 1.;
00377    fprintf(TecPlot,"%lf %lf %lf %lf %lf %lf\n",x,y,z,
00378       Plot[0][(vx*NBin+vy)*NBin+vz]*Norm0,Plot[1][(vx*NBin+vy)*NBin+vz]*Norm1,Plot[2][(vx*NBin+vy)*NBin+vz]*Norm2);
00379    // fprintf(TecPlot,"%lf %lf %lf %lf\n",x,y,z,Plot[0][vx*NBin+vy]*Norm0);
00380       }
00381     }
00382   }
00383   fclose(TecPlot);
00384   for(int t=0;t<NType;t++)
00385     free(Plot[t]);
00386   free(Plot);
00387   free(Count);
00388 }
00389 void ElPoly::HeaderPovRay(){
00390   FILE *HeaderPov = fopen("PovHeader.inc","w");
00391   double SphRad = 0.005;
00392   double CylRad = 0.002;
00393   double PerspAngle = 15;
00394   int ImSize[2] = {1000,1000};
00395   fprintf(HeaderPov,"// POV 3.x input script : plot.pov\n// command: povray +W%d +H%d -ISnap.pov -OSnap.pov.tga +P +X +A +FT +C\n",ImSize[0],ImSize[1]);
00396   fprintf(HeaderPov,"#if (version < 3.5)\n#error \"POV3DisplayDevice has been compiled for POV-Ray 3.5 or above.\\nPlease upgrade POV-Ray.\"\n#end\n");
00397   //include
00398     fprintf(HeaderPov,"// #include \"colors.inc\"\n   \
00399 //  #include \"stones.inc\"\n\
00400 //  #include \"textures.inc\"\n\
00401 //  #include \"shapes.inc\"\n\
00402 //  #include \"glass.inc\"\n\
00403 //  #include \"metals.inc\"\n\
00404 //  #include \"woods.inc\"\n");
00405   fprintf(HeaderPov,"\
00406 #declare clip_on=array[3] {0, 0, 0};\n\
00407 #declare clip=array[3];\n\
00408 #declare scaledclip=array[3];\n\
00409 #declare SphRad=%lf;\n\
00410 #declare CylRad=%lf;\n\
00411 #declare line_width=0.0020;\n",SphRad,CylRad);
00412   // //materials
00413   // fprintf(HeaderPov,"\
00414   // #declare Fotoni  = photons {\n\
00415   //       target\n\
00416   //       refraction on\n\
00417   //       reflection on\n\
00418   //   }\n");
00419   //define a point
00420   fprintf(HeaderPov,"\
00421 #macro ElPoint (P1, R1, C1)\n\
00422   #local T = texture { finish { ambient 1.0 diffuse 0.0 phong 0.0 specular 0.0 } pigment { C1 } }\n \
00423   #if(clip_on[2])\n\
00424   intersection {\n\
00425     sphere {P1, R1 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00426     clip[2]\n\
00427   }\n\
00428   #else\n\
00429   sphere {P1, R1 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00430   #end\n\
00431 #end\n");
00432   //define a line
00433   fprintf(HeaderPov,"\
00434 #macro ElLine (P1, P2, C1)\n\
00435   #local T = texture { finish { ambient 1.0 diffuse 0.0 phong 0.0 specular 0.0 } pigment { C1 } }\n\
00436   #if(clip_on[2])\n\
00437   intersection {\n\
00438     cylinder {P1, P2, line_width texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00439     clip[2]\n\
00440   }\n\
00441   #else\n\
00442   cylinder {P1, P2, line_width texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00443   #end\n\
00444 #end\n");
00445   //define a sphere
00446   fprintf(HeaderPov,"#macro ElSphere(P1, R1, C1)\n\
00447   #local T = texture { pigment { C1 } }\n\
00448   #local M = material{\n       \
00449     texture{\n\
00450       pigment{gradient y color_map{[0.4 C1][0.4 C1]}}\n\
00451       finish{ambient 0 diffuse 0.4 specular 1 roughness 0.0001 reflection 0.25}\n\
00452     }\n\
00453     interior{ior 1.33}\n            \
00454   }\n\
00455   #if(clip_on[2])\n\
00456   intersection {\n\
00457     sphere {P1, R1 material {M} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00458     clip[2]\n\
00459   }\n\
00460   #else\n\
00461   sphere {P1, R1 material {M} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00462   #end\n\
00463 #end\n");
00464   //define a cylinder
00465   fprintf(HeaderPov,"\
00466 #macro ElCylinder (P1, P2, R1, C1, O1)\n\
00467   #local T = texture { pigment { C1 } }\n\
00468   #if(clip_on[2])\n\
00469   intersection {\n\
00470     cylinder {P1, P2, R1 #if(O1) open #end texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow\n\
00471 }\n\
00472     clip[2]\n\
00473   }\n\
00474   #else\n\
00475   cylinder {P1, P2, R1 #if(O1) open #end texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00476   #end\n\
00477 #end\n");
00478   //define a triangle
00479   fprintf(HeaderPov,"\
00480 #macro ElTriangle (P1, N1, P2, N2, P3, N3, C1)\n\
00481   #local T = texture { pigment { C1 } }\n\
00482   smooth_triangle {P1, N1, P2, N2, P3, N3 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
00483 #end\n");
00484   //set the camera
00485   double Loc[3] = {0.5,4.0,4.};
00486   double LAt[3] = {0.5,0.5,0.5};
00487   double Up[3] = {0.,6.,0.};
00488   double Right[3] = {4.8,0.,0.};
00489   double Dir[3] = {0.,0.,4.};
00490   double Light0[3] = {0.5,4.0,4.0};
00491   double Light1[3] = {0.5,4.0,2.0};
00492   double CBack[3] = {1.,1.,1.};
00493   #ifdef __glut_h__
00494   // Draw *Dr;
00495   // Loc[0] = Dr->xp; Loc[1] = Dr->yp; Loc[2] = Dr->zp+Dr->zw;
00496   // Light0[0] = Dr->xl0; Light0[1] = Dr->yl0; Light0[2] = Dr->zl0;
00497   // Light1[0] = Dr->xl1; Light1[1] = Dr->yl1; Light1[2] = Dr->zl1;
00498   // CBack[0] = Dr->Rback; CBack[1] = Dr->Gback; CBack[2] = Dr->Bback;
00499   #endif //__glut_h__
00500   fprintf(HeaderPov,"camera {\n\
00501 //perspective orthographic fisheye ultra_wide_angle...\n\
00502   location <%.3f, %.3f, %.3f>\n\
00503   look_at <%.3f, %.3f, %.3f>\n\
00504   sky <0. 0. 1.>\n\
00505 //  up y\n\
00506 //  right 1.33*x\n\
00507 //  direction <%.3f, %.3f, %.3f>\n\
00508 //  translate <0.0, 1.0, 1.0>\n\
00509 //  rotate <0.0, 1.0, 1.0>\n\
00510     angle %f\n        \
00511 }\n",
00512      Loc[0],Loc[1],Loc[2],
00513      LAt[0],LAt[1],LAt[2],
00514      Dir[0],Dir[1],Dir[2],
00515      PerspAngle);
00516   //set the light0
00517   fprintf(HeaderPov,"\
00518 light_source {\n\
00519   <%lf, %lf, %lf>\n\
00520   color rgb<1.000, 1.000, 1.000>\n\
00521   parallel\n\
00522   point_at <.5, .5, .5>\n\
00523 }\n",Light0[0],Light0[1],Light0[2]);
00524   //set the light1
00525   fprintf(HeaderPov,"\
00526 light_source {\n\
00527   <%lf, %lf, %lf>\n\
00528   color rgb<1.000, 1.000, 1.000>\n\
00529   parallel\n\
00530   point_at <.5, .5, .5>\n\
00531 }\n",Light1[0],Light1[1],Light1[2]);
00532   //set the fog
00533   fprintf(HeaderPov,"\
00534   fog { fog_type   2\n\
00535       distance   50\n\
00536       color       rgb<1.000, 1.000, 1.000>\n\
00537       fog_offset 0.1\n\
00538       fog_alt    1.5\n\
00539       turbulence 1.8\n\
00540     }\n");
00541   //background
00542   fprintf(HeaderPov,"\
00543 background {\n\
00544   color rgb<%lf, %lf, %lf>\n\
00545 }\n",CBack[0],CBack[1],CBack[2]);
00546   //texture
00547   fprintf(HeaderPov,"\
00548 #default { texture {\n\
00549  finish { ambient 0.000 diffuse 0.650 phong 0.1 phong_size 40.000 specular 0.500 }\n\
00550 } }\n");
00551 }
00552 void ElPoly::ConvLattice(int NSample,char *FName){
00553   int NType = 3;
00554   double **Plot = (double **) calloc(3,sizeof(double));
00555   for(int t=0;t<NType;t++){
00556     Plot[t] = (double *)calloc(NSample*NSample,sizeof(double));
00557   }
00558   LoadDensFile(Plot,NSample);
00559   FILE *FWrite = fopen(FName,"w");
00560   fprintf(FWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(0),pEdge(1),pEdge(2),NSample,ChooseDraw(EL_PART));
00561   double InvNSample = 1./(double)NSample;
00562   for(int t=0;t<3;t++){
00563     for(int sx=0;sx<NSample;sx++){
00564       for(int sy=0;sy<NSample;sy++){
00565    //if(fabs(Plot[t][sx*NSample+sy]) < 0.01) continue;
00566    double x = sx*InvNSample*pEdge(CLat1);
00567    double y = sy*InvNSample*pEdge(CLat2);
00568    fprintf(FWrite,"{t[%d %d %d] x(%lf %lf %lf)}\n",t*pNPart()+sx*NSample+sy,t,t,x,y,Plot[t][sx*NSample+sy]);
00569       }
00570     }
00571   }
00572   for(int t=0;t<3;t++) free(Plot[t]);
00573   free(Plot);
00574 }