initwin.cpp

pic 模拟程序！面向对象

//  This initializes the windows for diagnostics.

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <txc_streams.h>
#include "sptlrgn.h"

#include <stdlib.h>
#include "oopiclist.h"
#include "ovector.h"
#include "diagn.h"
#include "ptclgrp.h"
#include "globalvars.h"
#include "boundary.h"
#include "gpdist.h"
#include "history.h"

extern "C" {
#include <math.h>
#include <stdio.h>
  //#include <malloc.h>

#ifdef NOX
#include <xgmini.h>
#else
#include <xgrafix.h>
void WinSet2D(char *PlotType, char *X_Label, char *Y_Label, char *W_Label,
             char *State, int ulx, int uly, Scalar X_Scale, Scalar Y_Scale,
             int X_Auto_Rescale, int Y_Auto_Rescale, Scalar X_Min,
             Scalar X_Max, Scalar Y_Min, Scalar Y_Max);
#endif
}

#define True 1
#define False 0

//  These are xgrafix flags!
//  They turn on/off the copying of the diagnostics from the pic
//  structures to the arrays that xgrafix must have.


void Diagnostics::InitWin(void) {
  // initialize the U ostrings 
  ostring U1, U2, U3;
  ostring U="u";
  //  We'll put the diagnostic calls to xgrafix here.

  jm = theSpace->getJ()+1;
  km = theSpace->getK()+1;
  km2 = km/2+1;

  jmNGD = theSpace->getJ();
  kmNGD = theSpace->getK();

  ostring prepend="";
  int i, isp;
	if(RegionName != (ostring)"Noname" && RegionName != (ostring)"") 
		prepend = RegionName + (ostring)": ";


  switch(theSpace->get_grid()->query_geometry())
	 {
	 case ZRGEOM:
		{
		  X1="z";
		  X2="r";
		  X3="phi";
		  K1="kz";
		  K2="kr";
                  U1=U+X1;
                  U2=U+X2;
                  U3=U+X3;
		  E1= prepend + "Ez";
		  E2= prepend + "Er";
		  E3= prepend + "Ephi";
		  B1= prepend + "Bz";
		  B2= prepend + "Br";
		  B3= prepend + "Bphi";
		  I1= prepend + "Iz";
		  I2= prepend + "Ir";
		  I3= prepend + "Iphi";
		  break;
		}
	 case ZXGEOM:
		{
		  X1="x";
		  X2="y";
		  X3="z";
		  K1="kx";
		  K2="ky";
                  U1=U+X1;
                  U2=U+X2;
                  U3=U+X3;
		  E1= prepend + "Ex";
		  E2= prepend + "Ey";
		  E3= prepend + "Ez";
		  B1= prepend + "Bx";
		  B2= prepend + "By";
		  B3= prepend + "Bz";
		  I1= prepend + "Ix";
		  I2= prepend + "Iy";
		  I3= prepend + "Iz";
		  break;
		}
	 };			
  char buf1[100];
  sprintf(buf1,"u%s vs %s, all species",X1.c_str(),X1.c_str());
  U1vsX1 =   prepend + buf1;
  sprintf(buf1,"u%s vs %s, all species",X2.c_str(),X1.c_str());
  U2vsX1 =   prepend + buf1;
  sprintf(buf1,"u%s vs %s, all species",X3.c_str(),X1.c_str());
  U3vsX1 =   prepend + buf1;
  sprintf(buf1,"u%s vs %s, all species",X1.c_str(),X2.c_str());
  U1vsX2=   prepend + buf1;
  sprintf(buf1,"u%s vs %s, all species",X2.c_str(),X2.c_str());
  U2vsX2=   prepend + buf1;
  sprintf(buf1,"u%s vs %s, all species",X3.c_str(),X2.c_str());
  U3vsX2=   prepend + buf1;

#ifndef NOX
  //  Set up the diagnostics 
  //----------------------------------------------------------------------//
  // Set up the new diagnostics

  DiagList *dlist = theSpace->getDiagList();
	oopicListIter<Diag> nextd(*dlist);
  for(nextd.restart();!nextd.Done(); nextd++) 
		nextd.current()->initwin();

  //----------------------------------------------------------------------//
  //set up electric fields
  XGSet3D( "linlinlin",X1,X2,E1,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)E, &jm, &km, 1 ,1,0,1,0,3,0);

  XGSet3D( "linlinlin",X1,X2,E2,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)E, &jm, &km, 1 ,1,0,1,0,3,1);

  XGSet3D( "linlinlin",X1,X2,E3,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)E, &jm, &km, 1 ,1,0,1,0,3,2);

#ifdef DEBUG_FIELDS
  XGSet3D( "linlinlin",X1,X2,"intEdl1",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intEdl, &jm, &km, 1 ,1,0,1,0,3,0);

  XGSet3D( "linlinlin",X1,X2,"intEdl2",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intEdl, &jm, &km, 1 ,1,0,1,0,3,1);

  XGSet3D( "linlinlin",X1,X2,"intEdl3",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intEdl, &jm, &km, 1 ,1,0,1,0,3,2);
#endif
  //set up magnetic fields

  XGSet3D( "linlinlin",X1,X2,B1,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)B, &jm, &km, 1,1,0,1,0,3,0 );

  XGSet3D( "linlinlin",X1,X2,B2,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)B, &jm, &km, 1,1,0,1,0,3,1 );

  XGSet3D( "linlinlin",X1,X2,B3,45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)B, &jm, &km, 1,1,0,1,0,3,2 );

#ifdef DEBUG_FIELDS
  XGSet3D( "linlinlin",X1,X2,"intBdS1",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intBdS, &jm, &km, 1 ,1,0,1,0,3,0);

  XGSet3D( "linlinlin",X1,X2,"intBdS2",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intBdS, &jm, &km, 1 ,1,0,1,0,3,1);

  XGSet3D( "linlinlin",X1,X2,"intBdS3",45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurfVector( x1_array,x2_array,(Scalar **)intBdS, &jm, &km, 1 ,1,0,1,0,3,2);
#endif

  //set up currents
  if(!electrostaticFlag) {
	 XGSet3D( "linlinlin",X1,X2,I1,45.0,225.0,"closed",1,1,
				1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
	 XGSurfVector( x1_array,x2_array,(Scalar **)I, &jm, &km, 1,1,0,1,0,3,0 );

	 XGSet3D( "linlinlin",X1,X2,I2,45.0,225.0,"closed",1,1,
				1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
	 XGSurfVector( x1_array,x2_array,(Scalar **)I, &jm, &km, 1,1,0,1,0,3,1 );

	 XGSet3D( "linlinlin",X1,X2,I3,45.0,225.0,"closed",1,1,
				1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
	 XGSurfVector( x1_array,x2_array,(Scalar **)I, &jm, &km, 1,1,0,1,0,3,2 );
  }
  //set up energies

  XGSet3D( "linlinlin",X1,X2,strdup((prepend +"Ue")),45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurf( x1_array,x2_array,Ue, &jm, &km, 1 );

  XGSet3D( "linlinlin",X1,X2,strdup((prepend +"Ub")),45.0,225.0,"closed",1,1,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurf( x1_array,x2_array,Ub, &jm, &km, 2 );

  // Set up poynting vector
  if(!electrostaticFlag)
	 {
		XGSetVec("vecvec",X1,X2,strdup(prepend+"Poynting Vector"),"closed",1,1,1.0
					,1.0,0,0,x1min,x1max,x2min,x2max);

		XGVectorVector(x1_array,x2_array,(Scalar **)S_array,
							(Scalar **)S_array,&jm,&km,2,1,0,1,0,3,0,3,1);

	 }


  //  set up rho
	char buffer[100];
	if (boltzmannFlag)
		sprintf(buffer,"PIC ");
	else 
		sprintf(buffer,"%s","");

  XGSet3D( "linlinlin",X1,X2,strdup(prepend + buffer + "rho"),45.0,250.0,"closed",1,400,
			 1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurf( x1_array,x2_array,rho, &jm, &km, 3 );

  /** 
   *  Plot the neutral gas densities, dad: 01/24/2001
   */ 
  oopicListIter<NGD> NGDIter(*ptrNGDList);
  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ) {
    XGSet3D( "linlinlin",X1,X2,
               strdup(prepend + "Neutral Gas Density for " + (NGDIter.current())->getGasType()),
               45.0,250.0,"closed",1,400,
               1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
    XGSurf( x1_arrayNGD,x2_arrayNGD, (NGDIter.current())->getNGDdata(), &jmNGD, &kmNGD, 3 );
  }
  
  //  set up boltzmann rho
  if (boltzmannFlag)
    {
      XGSet3D( "linlinlin",X1,X2,strdup(prepend +"Total Rho"),45.0,250.0,"closed",1,1,
               1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
      XGSurf( x1_array,x2_array,totalRho, &jm, &km, 3 );
    }
  if(electrostaticFlag) {
    //  set up phi
    XGSet3D( "linlinlin",X1,X2,strdup(prepend +"phi"),45.0,250.0,"closed",1,1,
             1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
    XGSurf( x1_array,x2_array,phi, &jm, &km, 3 );
  }
#ifdef DEBUG_PHI
  XGSet3D( "linlinlin",X1,X2,strdup(prepend +"phi error"),45.0,250.0,"closed",1,1,
           1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
  XGSurf( x1_array,x2_array,phi_err, &jm, &km, 3 );
#endif
  //  set up divderror
  if(!electrostaticFlag) {
    XGSet3D( "linlinlin",X1,X2,strdup(prepend +"Divergence Error"),45.0,250.0,"closed",1,1,
             1.0,1.0,1.0,0,0,1,x1min,x1max,x2min,x2max,0.0,1.0);
    XGSurf( x1_array,x2_array,divderror, &jm, &km, 3 );
  }
  
  //  The phase space diagnostics
  //	X2 vs X1
  for(isp=0;isp<number_of_species;isp++) {
    char * buf = new char[80];
    sprintf(buf,"%s-%s phase space for %s ",X1.c_str(),X2.c_str(),theSpecies[isp].name);
    
    WinSet2D( "linlin", X1, X2, strdup(prepend +buf),"open"
              ,700,400,1,1,0,0,x1min,x1max,x2min,x2max);
    delete [] buf; 
    
    XGScat2D(theSpecies[isp].x1positions,theSpecies[isp].x2positions,
             &(theSpecies[isp].nparticles_plot),isp);
    //draw any boundaries or structures
    oopicListIter<Boundary> nextb(*blist);
    for(nextb.restart();!nextb.Done(); nextb++) {
      Scalar x1=0.;
      Scalar x2=0.;
      Scalar y1=0.;
      Scalar y2=0.;
      int j1,j2,k1,k2;
      
      if (nextb.current()->queryFill()){
        Boundary ***NB;
        j1=nextb.current()->getMinJ();
        j2=nextb.current()->getMaxJ();
        k1=nextb.current()->getMinK();
        k2=nextb.current()->getMaxK();
        NB=nextb.current()->getFields()->get_grid()->GetCellMask();
        
        int j,k;
        int on=0;
        for (j=j1; j<j2; j++){
          for (k=k1; k<k2; k++){
            if (NB[j][k]==nextb.current()){
              if (!on){
                on=1;
                x1 = x1_array[j]; 
                y1 = x2_array[k];
              }
            }
            else 
              if (on) {
                on=0;
                int bcolor = nextb.current()->getBoundaryType();
                Scalar *fpoints = (Scalar *) calloc(32,sizeof(Scalar));
                x2 = x1_array[j+1];
                y2 = x2_array[k];
                fpoints[0]=x1;fpoints[1]=y1;fpoints[2]=x1;fpoints[3]=y2;
                fpoints[4]=x2;fpoints[5]=y2;fpoints[6]=x2;fpoints[7]=y1;
                fpoints[8]=x1;fpoints[9]=y1;
                XGStructureArray(5,FILLED,-bcolor,-bcolor,fpoints);
                free(fpoints);
              }
          }
          if (on){
            on = 0;
            Scalar *fpoints = (Scalar *) calloc(32,sizeof(Scalar));
            int bcolor = nextb.current()->getBoundaryType();
            fpoints[0]=x1;fpoints[1]=y1;fpoints[2]=x1;fpoints[3]=y2;
            fpoints[4]=x2;fpoints[5]=y2;fpoints[6]=x2;fpoints[7]=y1;
            fpoints[8]=x1;fpoints[9]=y1;
            x2 = x1_array[j+1];
            y2 = x2_array[k2];
            XGStructureArray(5,FILLED,-bcolor,-bcolor,fpoints);
            free(fpoints);
          }
        }
      }
      
      oopicListIter <LineSegment> lsI(*nextb.current()->getSegmentList());
      for(lsI.restart();!lsI.Done();lsI++) {
        int j1,j2,k1,k2,normal,*points;  //local copies of above
        j1=(int)lsI.current()->A.e1();
        k1=(int)lsI.current()->A.e2();
        j2=(int)lsI.current()->B.e1();
        k2=(int)lsI.current()->B.e2();
        points=lsI.current()->points;
        normal=lsI.current()->normal;
        
        x1 = x1_array[j1]; 
        x2 = x1_array[j2];
        y1 = x2_array[k1];
        y2 = x2_array[k2];
        int bcolor = nextb.current()->getBoundaryType();
        
        // Different sorts of boundaries need to be drawn differently
        switch(bcolor) {
        case DIELECTRIC_TRIANGLE_UP: // positive-normal DielectricTriangle
          {				 
            Scalar * fpoints = (Scalar *) calloc(4*abs(j2-j1+3),sizeof(Scalar));
            int j;
            for(j=0;j<4*(j2-j1)+4;j+=2)
              { fpoints[j] = x1_array[points[j]];
              fpoints[j+1] = x2_array[points[j+1]];
              }
            fpoints[j] = (y1 < y2)? x1 :x2;
            fpoints[j+1] = MAX(y1,y2);
            fpoints[j+2]= fpoints[0];
            fpoints[j+3]= fpoints[1];
            XGStructureArray(2*(j2-j1)+4,FILLED,bcolor - 100, bcolor -100,
                             fpoints);
            
            //				delete [] points;
            free(fpoints);
            
            break;
          }
        case DIELECTRIC_TRIANGLE_DOWN: //negative-normal DielectricTriangle
          {
            Scalar * fpoints = (Scalar *) calloc(4*abs(j2-j1+3),sizeof(Scalar));
            int *points = segmentate(j1,k1,j2,k2);
            int j;
            for(j=0;j<4*(j2-j1)+4;j+=2)