C++ ComplexMatrix::Rlo方法代码示例

ComplexMatrix MatrixSqrtHermitean (const ComplexMatrix& A)
//
// Returns the matrix square root of a hermitean matrix A.
// The argument matrix is not checked for hermiticity !
//
//
{  
    // get the dimension information
    int lo = A.Clo(),
        hi = A.Chi();

    // columns and rows must have the same range
    if (A.Rlo() != lo || A.Rhi() != hi) 
        Matpack.Error(Mat::UnspecifiedError,"MatrixSqrtHermitean: hermitean matrix must be square"); 
    
    Matrix z(lo,hi,lo,hi), 
           zr(lo,hi,lo,hi),
           zi(lo,hi,lo,hi);

    Vector d(lo,hi);

    int i,j;

    // create packed form of hermitean matrix A in matrix z
    for (i = lo; i <= hi; i++) {        
        for (j = lo; j < i; j++) {
            z[j][i]=imag(A[i][j]);
            z[i][j]=real(A[i][j]);
        }
        z[i][i]=real(A[i][i]);
    }
    
    // diagonalize z
    EigenSystemHermitean(z,d,zr,zi,false,30);

    // combine real and imaginary parts of eigenvectors (zr,zi) in complex matrix
    ComplexMatrix U(zr,zi);
  
    ComplexVector dc(lo,hi);

    for (i = lo; i <= hi; i++) {        
        if (d[i] < 0) 
          dc[i]=complex<double>(0,sqrt(-d[i]));
        else 
          dc[i]=complex<double>(sqrt(d[i]),0);
    }

    // calculate matrix square root
    U = U * Diagonal(dc) * U.Hermitean();

    return U.Value();    
}

int main (void)
{
  MpImage image(size,size);		     // allocate size x size raster image

  MpImage backimg;			                 // read background image
  if (back_image) {
    if ( ! backimg.ReadAnyFile(back_image) )
      Matpack.Error("Can't read \"%s\"",back_image);
  }

  // read potential matrix
  ifstream pot_stream("pot");
  if ( ! pot_stream) Matpack.Error("Can't open pot data file \"%s\"","pot");
  pot_stream >> pot;
  pot_stream.close();

  for (int i = 0; i <= 175; i += 5) {  			   // loop through frames

    fprintf(stdout,"\rframe %d...",i); fflush(stdout);         // what's going on

    // read complex wave function
    char file_name[200], pipe_name[200], output_name[200];
    sprintf(file_name,"psi%d.gz",i);		           // generate input data file name
    sprintf(pipe_name,"gunzip -c %s > psi.dat",file_name); // decompress data file 
    sprintf(output_name,"psi%03d.jpg",i);                  // generate output file name

    system(pipe_name);
    ifstream data_stream("psi.dat");
    if ( ! data_stream) Matpack.Error("Can't open data file psi.dat");
    data_stream >> psi;
    data_stream.close();
    unlink("psi.dat");

    // consistency checks
    if ( pot.Rlo() != psi.Rlo() || pot.Rhi() != psi.Rhi() || 
	 pot.Clo() != psi.Clo() || pot.Chi() != psi.Chi() )
      Matpack.Error("non-conformant index range of potential and wave function");

    Scene scene;				   // define scene for 3d drawing

    // create surface
    double u0 = psi.Clo(), u1 = psi.Chi(), 
           v0 = psi.Rlo(), v1 = psi.Rhi();
    int    nu = psi.Chi()-psi.Clo(), nv = psi.Rhi()-psi.Rlo(),  
           su = 0, sv = 0,  periodic = 0;
    ParametricSurface(scene, height_fcn, u0,u1,nu,su, v0,v1,nv,sv, 
		      periodic, Identity, color_fcn);

    scene.BoxRatios(1,z_aspect,1);               // scale scene to box with x:y:z

    float dist = 1.6;			                    // distance parameter
    scene.Look(Vector3D( 1.1*dist, 1.5*dist, 1.4*dist), // camera position vector
	       Vector3D(-1.1*dist,-1.65*dist,-1.4*dist), // look direction vector
	       FieldOfView(45), 		 // field of view angle in degree
	       0);		             // "twist your head" angle in degree

    scene.SetGlobalShading(Facet::Gouraud);   
                                       // shading (None, Flat, Gouraud, Phong...)
    scene.SetHiddenSurface(preview ? Scene::DepthSort : Scene::TopoSort);

    // edgeline drawing option (None,...)
    scene.SetGlobalEdgeLines(show_grid ? Facet::Individual : None);

    // Setting per-vertex coloring is the best choice to get smoothly changing
    // colors on the surface. This requires Gouraud or Phong shading to be
    // switched on. If you use per-facet coloring then the facet edges are still
    // slightly visible ("Scene::PerFacet")
    scene.SetColoring(Scene::PerVertex);

    // Define ambient light, RGB color, and specular light and exponent
    // If per-vertex coloring or per-facet coloring is set, then only the ambient
    // light factor and the specular light factor and exponent are in effect, and
    // the diffuse color is ignored (the vertex or facet color is used instead)
    Material material(0.6, ColorF(1.,0.8,0.), 0.3,1.8); 
    scene.SetFacetMaterial(material);		    // set all facets to material

    scene.Open(image);			        // direct drawing to raster image
    if (back_image)
      image.InsertTiled(backimg);
    else
      scene.SetBackground(back_color); 	        // draw background with RGB color
    scene.SetColor(edge_color);		           // define color for edge lines
    scene.Show();			  		    // render the surface
    scene.Close();				       // call always after close

    image.WriteJpegFile(output_name,jpeg_quality,jpeg_smooth);// write JPEG image

    // write GIF image, if you prefer that
    //image.WriteGifFile(output_name);     // write image as GIF

  }

  cout << "\nfinished" << endl;
}