C++ Parameters类代码示例

FuncDeclaration *StructDeclaration::buildXopEquals(Scope *sc)
{
    if (!search_function(this, Id::eq))
        return NULL;

    /* static bool__xopEquals(in void* p, in void* q) {
     *     return ( *cast(const S*)(p) ).opEquals( *cast(const S*)(q) );
     * }
     */

    Parameters *parameters = new Parameters;
    parameters->push(new Parameter(STCin, Type::tvoidptr, Id::p, NULL));
    parameters->push(new Parameter(STCin, Type::tvoidptr, Id::q, NULL));
    TypeFunction *tf = new TypeFunction(parameters, Type::tbool, 0, LINKd);
    tf = (TypeFunction *)tf->semantic(loc, sc);

    Identifier *id = Lexer::idPool("__xopEquals");
    FuncDeclaration *fop = new FuncDeclaration(loc, 0, id, STCstatic, tf);

    Expression *e = new CallExp(0,
        new DotIdExp(0,
            new PtrExp(0, new CastExp(0,
                new IdentifierExp(0, Id::p), type->pointerTo()->constOf())),
            Id::eq),
        new PtrExp(0, new CastExp(0,
            new IdentifierExp(0, Id::q), type->pointerTo()->constOf())));

    fop->fbody = new ReturnStatement(loc, e);

    size_t index = members->dim;
    members->push(fop);

    sc = sc->push();
    sc->stc = 0;
    sc->linkage = LINKd;

    unsigned errors = global.startGagging();
    fop->semantic(sc);
    if (errors == global.gaggedErrors)
    {   fop->semantic2(sc);
        if (errors == global.gaggedErrors)
        {   fop->semantic3(sc);
            if (errors == global.gaggedErrors)
                fop->addMember(sc, this, 1);
        }
    }
    if (global.endGagging(errors))    // if errors happened
    {
        members->remove(index);

        if (!xerreq)
        {
            Expression *e = new IdentifierExp(loc, Id::empty);
            e = new DotIdExp(loc, e, Id::object);
            e = new DotIdExp(loc, e, Lexer::idPool("_xopEquals"));
            e = e->semantic(sc);
            Dsymbol *s = getDsymbol(e);
            FuncDeclaration *fd = s->isFuncDeclaration();

            xerreq = fd;
        }
        fop = xerreq;
    }

    sc->pop();

    return fop;
}

  ToneProfile::ToneProfile(tone_profile_t whichProfile, bool majorScale, const Parameters& params){

    float p[12];
    if(whichProfile == TONE_PROFILE_SILENT){
      p[0]=0; p[1]=0;
      p[2]=0; p[3]=0;
      p[4]=0; p[5]=0;
      p[6]=0; p[7]=0;
      p[8]=0; p[9]=0;
      p[10]=0; p[11]=0;
    }else if(whichProfile == TONE_PROFILE_TEMPERLEY){
      if(majorScale){
        p[0]=5.0; p[1]=2.0;
        p[2]=3.5; p[3]=2.0;
        p[4]=4.5;
        p[5]=4.0; p[6]=2.0;
        p[7]=4.5; p[8]=2.0;
        p[9]=3.5; p[10]=1.5;
        p[11]=4.0;
      }else{
        p[0]=5.0; p[1]=2.0;
        p[2]=3.5;
        p[3]=4.5; p[4]=2.0;
        p[5]=4.0; p[6]=2.0;
        p[7]=4.5;
        p[8]=3.5; p[9]=2.0;
        p[10]=1.5; p[11]=4.0;
      }
    }else if(whichProfile == TONE_PROFILE_GOMEZ){
      if(majorScale){
        p[0]=0.82; p[1]=0.00;
        p[2]=0.55; p[3]=0.00;
        p[4]=0.53;
        p[5]=0.30; p[6]=0.08;
        p[7]=1.00; p[8]=0.00;
        p[9]=0.38; p[10]=0.00;
        p[11]=0.47;
      }else{
        p[0]=0.81; p[1]=0.00;
        p[2]=0.53;
        p[3]=0.54; p[4]=0.00;
        p[5]=0.27; p[6]=0.07;
        p[7]=1.00;
        p[8]=0.27; p[9]=0.07;
        p[10]=0.10; p[11]=0.36;
      }
    }else if(whichProfile == TONE_PROFILE_SHAATH){
      if(majorScale){
        p[0]=6.6; p[1]=2.0;
        p[2]=3.5; p[3]=2.3;
        p[4]=4.6;
        p[5]=4.0; p[6]=2.5;
        p[7]=5.2; p[8]=2.4;
        p[9]=3.7; p[10]=2.3;
        p[11]=3.4;
      }else{
        p[0]=6.5; p[1]=2.7;
        p[2]=3.5;
        p[3]=5.4; p[4]=2.6;
        p[5]=3.5; p[6]=2.5;
        p[7]=5.2;
        p[8]=4.0; p[9]=2.7;
        p[10]=4.3; p[11]=3.2;
      }
    }else if(whichProfile == TONE_PROFILE_KRUMHANSL){
      if(majorScale){
        p[0]=6.35; p[1]=2.23;
        p[2]=3.48; p[3]=2.33;
        p[4]=4.38;
        p[5]=4.09; p[6]=2.52;
        p[7]=5.19; p[8]=2.39;
        p[9]=3.66; p[10]=2.29;
        p[11]=2.88;
      }else{
        p[0]=6.33; p[1]=2.68;
        p[2]=3.52;
        p[3]=5.38; p[4]=2.60;
        p[5]=3.53; p[6]=2.54;
        p[7]=4.75;
        p[8]=3.98; p[9]=2.69;
        p[10]=3.34; p[11]=3.17;
      }
    }else{ // Custom
      std::vector<float> ctp = params.getCustomToneProfile();
      if(majorScale){
        for (unsigned int i=0; i<12; i++)
          p[i] = (float)ctp[i];
      }else{
        for (unsigned int i=0; i<12; i++)
          p[i] = (float)ctp[i+12];
      }
    }

    // copy into doubly-linked circular list
    tonic = new Binode<float>(p[0]);
    Binode<float> *q = tonic;
    for (unsigned int i=1; i<12; i++){
      q->r = new Binode<float>(p[i]);
      q->r->l = q;
      q = q->r;
//.........这里部分代码省略.........

RadCalcer::RadCalcer( Parameters& params_, long opacity_, bool correlated_ )
  : zDist( params_, opacity_ ), qperps( opacity_, correlated_ ), n( opacity_ )
{
  std::vector<double> ReturnedParamsDouble;
  std::vector<long> ReturnedParamsLong;
  std::list<std::string> ReturnedParamsString;
  std::list<std::string>::iterator it;

  // First, we get the parameters of the medium
  // This is: 1st = mu, 2nd = temperature, 3rd = gluon mass, 4th = gluon mean free path
  // Here, we only need the gluon mass
  // the zDist will need the others
  ReturnedParamsDouble = params_.GetParametersDouble( "@mediumParams" );
  _mg = ReturnedParamsDouble[2];

  // Next, get the jet flavour
  ReturnedParamsString = params_.GetParametersString( "@jetFlavour" );
  std::string jetFlavour = ReturnedParamsString.front();
  if ( jetFlavour == "Gluon" )
  {
    _cr = 3.; _mass = _mg;
  }
  else if ( jetFlavour == "Light" )
  {
    _cr = 4./3.; _mass = _mg / sqrt(2);
  }
  else if ( jetFlavour == "Charm" )
  {
    _cr = 4./3.; _mass = 1.2;
  }
  else if ( jetFlavour == "Bottom" )
  {
    _cr = 4./3.; _mass = 4.75;
  }
  else
  {
    std::cerr << "@jetFlavour not understood" << std::endl;
    exit(0);
  }

  // Now, check whether also specifying a jet mass
  ReturnedParamsDouble = params_.GetParametersDouble( "@jetMassDirect" );
  if ( ReturnedParamsDouble.size() > 0 )
    _mass = ReturnedParamsDouble[0];
  
  // Now, the momentum, to give the jet energy
  ReturnedParamsDouble = params_.GetParametersDouble( "@jetMomentum" );
  double jetMomentum = ReturnedParamsDouble[0];
  // Now calculate and set energy, and mass
  _en = sqrt( _mass*_mass + jetMomentum*jetMomentum );

  // The jet path length in the medium
  ReturnedParamsDouble = params_.GetParametersDouble( "@pathLength" );
  _length = ReturnedParamsDouble[0];

  // The setting on the k max
  ReturnedParamsLong = params_.GetParametersLong( "@limitSet" );
  _switchkmax = ReturnedParamsLong[0];

  // The strong coupling, alpha_s
  ReturnedParamsString = params_.GetParametersString( "@alpha" );
  it = ReturnedParamsString.begin();
  if ( *it == "fixed" )
  {
    ++it;
    _alphas = boost::lexical_cast<double>( *it );
  }
  else
  {
    std::cerr << "Radcalcer, @alpha not understood";
    std::cerr << std::endl;
    exit(0);
  }

  ReturnedParamsString = params_.GetParametersString( "@incClassicalDiffusion" );
  if ( ReturnedParamsString.front() == "yes" )
    _diffexclude = false;
  else if ( ReturnedParamsString.front() == "no" )
    _diffexclude = true;
  else
  {
    std::cerr << "Radcalcer, @incClassicalDiffusion not understood as yes or no";
    std::cerr << std::endl;
    exit(0);
  }


  _correlated = correlated_;
}

FuncDeclaration *buildXopCmp(StructDeclaration *sd, Scope *sc)
{
    //printf("StructDeclaration::buildXopCmp() %s\n", toChars());
    if (Dsymbol *cmp = search_function(sd, Id::cmp))
    {
        if (FuncDeclaration *fd = cmp->isFuncDeclaration())
        {
            TypeFunction *tfcmpptr;
            {
                Scope scx;

                /* const int opCmp(ref const S s);
                 */
                Parameters *parameters = new Parameters;
                parameters->push(new Parameter(STCref | STCconst, sd->type, NULL, NULL));
                tfcmpptr = new TypeFunction(parameters, Type::tint32, 0, LINKd);
                tfcmpptr->mod = MODconst;
                tfcmpptr = (TypeFunction *)tfcmpptr->semantic(Loc(), &scx);
            }
            fd = fd->overloadExactMatch(tfcmpptr);
            if (fd)
                return fd;
        }
    }
    else
    {
#if 0   // FIXME: doesn't work for recursive alias this
        /* Check opCmp member exists.
         * Consider 'alias this', but except opDispatch.
         */
        Expression *e = new DsymbolExp(sd->loc, sd);
        e = new DotIdExp(sd->loc, e, Id::cmp);
        Scope *sc2 = sc->push();
        e = e->trySemantic(sc2);
        sc2->pop();
        if (e)
        {
            Dsymbol *s = NULL;
            switch (e->op)
            {
                case TOKoverloadset:    s = ((OverExp *)e)->vars;       break;
                case TOKimport:         s = ((ScopeExp *)e)->sds;       break;
                case TOKvar:            s = ((VarExp *)e)->var;         break;
                default:                break;
            }
            if (!s || s->ident != Id::cmp)
                e = NULL;   // there's no valid member 'opCmp'
        }
        if (!e)
            return NULL;    // bitwise comparison would work
        /* Essentially, a struct which does not define opCmp is not comparable.
         * At this time, typeid(S).compare might be correct that throwing "not implement" Error.
         * But implementing it would break existing code, such as:
         *
         * struct S { int value; }  // no opCmp
         * int[S] aa;   // Currently AA key uses bitwise comparison
         *              // (It's default behavior of TypeInfo_Strust.compare).
         *
         * Not sure we should fix this inconsistency, so just keep current behavior.
         */
#else
        return NULL;
#endif
    }

    if (!sd->xerrcmp)
    {
        // object._xopCmp
        Identifier *id = Lexer::idPool("_xopCmp");
        Expression *e = new IdentifierExp(sd->loc, Id::empty);
        e = new DotIdExp(sd->loc, e, Id::object);
        e = new DotIdExp(sd->loc, e, id);
        e = e->semantic(sc);
        Dsymbol *s = getDsymbol(e);
        if (!s)
        {
            ::error(Loc(), "ICE: %s not found in object module. You must update druntime", id->toChars());
            fatal();
        }
        assert(s);
        sd->xerrcmp = s->isFuncDeclaration();
    }

    Loc declLoc = Loc();    // loc is unnecessary so __xopCmp is never called directly
    Loc loc = Loc();        // loc is unnecessary so errors are gagged

    Parameters *parameters = new Parameters;
    parameters->push(new Parameter(STCref | STCconst, sd->type, Id::p, NULL));
    parameters->push(new Parameter(STCref | STCconst, sd->type, Id::q, NULL));
    TypeFunction *tf = new TypeFunction(parameters, Type::tint32, 0, LINKd);

    Identifier *id = Id::xopCmp;
    FuncDeclaration *fop = new FuncDeclaration(declLoc, Loc(), id, STCstatic, tf);

    Expression *e1 = new IdentifierExp(loc, Id::p);
    Expression *e2 = new IdentifierExp(loc, Id::q);
    Expression *e = new CallExp(loc, new DotIdExp(loc, e2, Id::cmp), e1);

    fop->fbody = new ReturnStatement(loc, e);

//.........这里部分代码省略.........

//.........这里部分代码省略.........
        VarDeclaration *vd = s->isVarDeclaration();
        if (vd && !vd->isDataseg())
        {
            if (vd->init)
            {
                // Should examine init to see if it is really all 0's
                zeroInit = 0;
                break;
            }
            else
            {
                if (!vd->type->isZeroInit(loc))
                {
                    zeroInit = 0;
                    break;
                }
            }
        }
    }

#if DMDV1
    /* This doesn't work for DMDV2 because (ref S) and (S) parameter
     * lists will overload the same.
     */
    /* The TypeInfo_Struct is expecting an opEquals and opCmp with
     * a parameter that is a pointer to the struct. But if there
     * isn't one, but is an opEquals or opCmp with a value, write
     * another that is a shell around the value:
     *  int opCmp(struct *p) { return opCmp(*p); }
     */

    TypeFunction *tfeqptr;
    {
        Parameters *arguments = new Parameters;
        Parameter *arg = new Parameter(STCin, handle, Id::p, NULL);

        arguments->push(arg);
        tfeqptr = new TypeFunction(arguments, Type::tint32, 0, LINKd);
        tfeqptr = (TypeFunction *)tfeqptr->semantic(Loc(), sc);
    }

    TypeFunction *tfeq;
    {
        Parameters *arguments = new Parameters;
        Parameter *arg = new Parameter(STCin, type, NULL, NULL);

        arguments->push(arg);
        tfeq = new TypeFunction(arguments, Type::tint32, 0, LINKd);
        tfeq = (TypeFunction *)tfeq->semantic(Loc(), sc);
    }

    Identifier *id = Id::eq;
    for (int i = 0; i < 2; i++)
    {
        Dsymbol *s = search_function(this, id);
        FuncDeclaration *fdx = s ? s->isFuncDeclaration() : NULL;
        if (fdx)
        {   FuncDeclaration *fd = fdx->overloadExactMatch(tfeqptr);
            if (!fd)
            {   fd = fdx->overloadExactMatch(tfeq);
                if (fd)
                {   // Create the thunk, fdptr
                    FuncDeclaration *fdptr = new FuncDeclaration(loc, loc, fdx->ident, STCundefined, tfeqptr);
                    Expression *e = new IdentifierExp(loc, Id::p);
                    e = new PtrExp(loc, e);
                    Expressions *args = new Expressions();

int main(int argInN, char* argIn[]) {
   
    time(&g_statsAll.timeStart);
   
    Parameters *P = new Parameters; //all parameters
       
    P->inputParameters(argInN, argIn);
    
    *(P->inOut->logStdOut) << timeMonthDayTime(g_statsAll.timeStart) << " ..... Started STAR run\n" <<flush;           
    
    //generate genome
    if (P->runMode=="genomeGenerate") {
        genomeGenerate(P);
        (void) sysRemoveDir (P->outFileTmp);        
        P->inOut->logMain << "DONE: Genome generation, EXITING\n" << flush;
        exit(0);
    } else if (P->runMode!="alignReads") {
        P->inOut->logMain << "EXITING because of INPUT ERROR: unknown value of input parameter runMode=" <<P->runMode<<endl<<flush;
        exit(1);
    };
    
    Genome mainGenome (P);
    mainGenome.genomeLoad();
    if (P->genomeLoad=="LoadAndExit" || P->genomeLoad=="Remove") 
    {
        return 0;
    };
    
    P->twoPass.pass2=false; //this is the 1st pass    
    
    SjdbClass sjdbLoci;

    if (P->sjdbInsert.pass1) 
    {
        Parameters *P1=new Parameters;
        *P1=*P;
        sjdbInsertJunctions(P, P1, mainGenome, sjdbLoci);
    };

    //calculate genome-related parameters
    Transcriptome *mainTranscriptome=NULL;
    
    
/////////////////////////////////////////////////////////////////////////////////////////////////START
    if (P->runThreadN>1) {
        g_threadChunks.threadArray=new pthread_t[P->runThreadN];
        pthread_mutex_init(&g_threadChunks.mutexInRead, NULL);
        pthread_mutex_init(&g_threadChunks.mutexOutSAM, NULL);
        pthread_mutex_init(&g_threadChunks.mutexOutBAM1, NULL);
        pthread_mutex_init(&g_threadChunks.mutexOutUnmappedFastx, NULL);
        pthread_mutex_init(&g_threadChunks.mutexOutFilterBySJout, NULL);
        pthread_mutex_init(&g_threadChunks.mutexStats, NULL);
        pthread_mutex_init(&g_threadChunks.mutexBAMsortBins, NULL);
    };

    g_statsAll.progressReportHeader(P->inOut->logProgress);    
    
    if (P->twoPass.yes) {//2-pass
        //re-define P for the pass1
        
        Parameters *P1=new Parameters;
        *P1=*P;
        //turn off unnecessary calculations
        P1->outSAMtype[0]="None";
        P1->outSAMbool=false;
        P1->outBAMunsorted=false;
        P1->outBAMcoord=false;
    
        P1->chimSegmentMin=0;
        
        P1->quant.yes=false;
        P1->quant.trSAM.yes=false;
        P1->quant.geCount.yes=false;
        
        P1->outFilterBySJoutStage=0;
        
        P1->outReadsUnmapped="None";
        
        P1->outFileNamePrefix=P->twoPass.dir;

        P1->readMapNumber=P->twoPass.pass1readsN;
//         P1->inOut->logMain.open((P1->outFileNamePrefix + "Log.out").c_str());

        g_statsAll.resetN();
        time(&g_statsAll.timeStartMap);
        P->inOut->logProgress << timeMonthDayTime(g_statsAll.timeStartMap) <<"\tStarted 1st pass mapping\n" <<flush;
        *P->inOut->logStdOut << timeMonthDayTime(g_statsAll.timeStartMap) << " ..... Started 1st pass mapping\n" <<flush;

        //run mapping for Pass1
        ReadAlignChunk *RAchunk1[P->runThreadN];        
        for (int ii=0;ii<P1->runThreadN;ii++) {
            RAchunk1[ii]=new ReadAlignChunk(P1, mainGenome, mainTranscriptome, ii);
        };    
        mapThreadsSpawn(P1, RAchunk1);
        outputSJ(RAchunk1,P1); //collapse and output junctions
//         for (int ii=0;ii<P1->runThreadN;ii++) {
//             delete [] RAchunk[ii];
//         };          
        
        time_t rawtime; time (&rawtime);
//.........这里部分代码省略.........

FuncDeclaration *buildOpAssign(StructDeclaration *sd, Scope *sc)
{
    if (FuncDeclaration *f = hasIdentityOpAssign(sd, sc))
    {
        sd->hasIdentityAssign = true;
        return f;
    }
    // Even if non-identity opAssign is defined, built-in identity opAssign
    // will be defined.

    if (!needOpAssign(sd))
        return NULL;

    //printf("StructDeclaration::buildOpAssign() %s\n", toChars());
    StorageClass stc = STCsafe | STCnothrow | STCpure | STCnogc;
    Loc declLoc = sd->loc;
    Loc loc = Loc();    // internal code should have no loc to prevent coverage

    if (sd->dtor || sd->postblit)
    {
        if (!sd->type->isAssignable())  // Bugzilla 13044
            return NULL;
        if (sd->dtor)
        {
            stc = mergeFuncAttrs(stc, sd->dtor);
            if (stc & STCsafe)
                stc = (stc & ~STCsafe) | STCtrusted;
        }
    }
    else
    {
        for (size_t i = 0; i < sd->fields.dim; i++)
        {
            VarDeclaration *v = sd->fields[i];
            if (v->storage_class & STCref)
                continue;
            Type *tv = v->type->baseElemOf();
            if (tv->ty == Tstruct)
            {
                TypeStruct *ts = (TypeStruct *)tv;
                if (FuncDeclaration *f = hasIdentityOpAssign(ts->sym, sc))
                    stc = mergeFuncAttrs(stc, f);
            }
        }
    }

    Parameters *fparams = new Parameters;
    fparams->push(new Parameter(STCnodtor, sd->type, Id::p, NULL));
    Type *tf = new TypeFunction(fparams, sd->handleType(), 0, LINKd, stc | STCref);

    FuncDeclaration *fop = new FuncDeclaration(declLoc, Loc(), Id::assign, stc, tf);

    Expression *e = NULL;
    if (stc & STCdisable)
    {
    }
    else if (sd->dtor || sd->postblit)
    {
        /* Do swap this and rhs
         *    tmp = this; this = s; tmp.dtor();
         */
        //printf("\tswap copy\n");
        Identifier *idtmp = Lexer::uniqueId("__tmp");
        VarDeclaration *tmp = NULL;
        AssignExp *ec = NULL;
        if (sd->dtor)
        {
            tmp = new VarDeclaration(loc, sd->type, idtmp, new VoidInitializer(loc));
            tmp->noscope = 1;
            tmp->storage_class |= STCtemp | STCctfe;
            e = new DeclarationExp(loc, tmp);
            ec = new BlitExp(loc, new VarExp(loc, tmp), new ThisExp(loc));
            e = Expression::combine(e, ec);
        }
        ec = new BlitExp(loc, new ThisExp(loc), new IdentifierExp(loc, Id::p));
        e = Expression::combine(e, ec);
        if (sd->dtor)
        {
            /* Instead of running the destructor on s, run it
             * on tmp. This avoids needing to copy tmp back in to s.
             */
            Expression *ec2 = new DotVarExp(loc, new VarExp(loc, tmp), sd->dtor, 0);
            ec2 = new CallExp(loc, ec2);
            e = Expression::combine(e, ec2);
        }
    }
    else
    {
        /* Do memberwise copy
         */
        //printf("\tmemberwise copy\n");
        for (size_t i = 0; i < sd->fields.dim; i++)
        {
            VarDeclaration *v = sd->fields[i];
            // this.v = s.v;
            AssignExp *ec = new AssignExp(loc,
                new DotVarExp(loc, new ThisExp(loc), v, 0),
                new DotVarExp(loc, new IdentifierExp(loc, Id::p), v, 0));
            e = Expression::combine(e, ec);
        }
//.........这里部分代码省略.........

FuncDeclaration *StructDeclaration::buildOpAssign(Scope *sc)
{
    if (!needOpAssign())
        return NULL;

    //printf("StructDeclaration::buildOpAssign() %s\n", toChars());

    FuncDeclaration *fop = NULL;

    Parameter *param = new Parameter(STCnodtor, type, Id::p, NULL);
    Parameters *fparams = new Parameters;
    fparams->push(param);
    Type *ftype = new TypeFunction(fparams, handle, FALSE, LINKd);
#if STRUCTTHISREF
    ((TypeFunction *)ftype)->isref = 1;
#endif

    fop = new FuncDeclaration(loc, 0, Id::assign, STCundefined, ftype);

    Expression *e = NULL;
    if (postblit)
    {   /* Swap:
         *    tmp = *this; *this = s; tmp.dtor();
         */
        //printf("\tswap copy\n");
        Identifier *idtmp = Lexer::uniqueId("__tmp");
        VarDeclaration *tmp;
        AssignExp *ec = NULL;
        if (dtor)
        {
            tmp = new VarDeclaration(0, type, idtmp, new VoidInitializer(0));
            tmp->noscope = 1;
            tmp->storage_class |= STCctfe;
            e = new DeclarationExp(0, tmp);
            ec = new AssignExp(0,
                new VarExp(0, tmp),
#if STRUCTTHISREF
                new ThisExp(0)
#else
                new PtrExp(0, new ThisExp(0))
#endif
                );
            ec->op = TOKblit;
            e = Expression::combine(e, ec);
        }
        ec = new AssignExp(0,
#if STRUCTTHISREF
                new ThisExp(0),
#else
                new PtrExp(0, new ThisExp(0)),
#endif
                new IdentifierExp(0, Id::p));
        ec->op = TOKblit;
        e = Expression::combine(e, ec);
        if (dtor)
        {
            /* Instead of running the destructor on s, run it
             * on tmp. This avoids needing to copy tmp back in to s.
             */
            Expression *ec = new DotVarExp(0, new VarExp(0, tmp), dtor, 0);
            ec = new CallExp(0, ec);
            e = Expression::combine(e, ec);
        }
    }
    else
    {   /* Do memberwise copy
         */
        //printf("\tmemberwise copy\n");
        for (size_t i = 0; i < fields.dim; i++)
        {
            Dsymbol *s = fields.tdata()[i];
            VarDeclaration *v = s->isVarDeclaration();
            assert(v && v->storage_class & STCfield);
            // this.v = s.v;
            AssignExp *ec = new AssignExp(0,
                new DotVarExp(0, new ThisExp(0), v, 0),
                new DotVarExp(0, new IdentifierExp(0, Id::p), v, 0));
            ec->op = TOKblit;
            e = Expression::combine(e, ec);
        }
    }
    Statement *s1 = new ExpStatement(0, e);

    /* Add:
     *   return this;
     */
    e = new ThisExp(0);
    Statement *s2 = new ReturnStatement(0, e);

    fop->fbody = new CompoundStatement(0, s1, s2);

    members->push(fop);
    fop->addMember(sc, this, 1);

    sc = sc->push();
    sc->stc = 0;
    sc->linkage = LINKd;

    fop->semantic(sc);

//.........这里部分代码省略.........

FuncDeclaration *StructDeclaration::buildOpEquals(Scope *sc)
{
    if (!needOpEquals())
        return NULL;
    //printf("StructDeclaration::buildOpEquals() %s\n", toChars());
    Loc loc = this->loc;

    Parameters *parameters = new Parameters;
#if STRUCTTHISREF
    // bool opEquals(ref const T) const;
    Parameter *param = new Parameter(STCref, type->constOf(), Id::p, NULL);
#else
    // bool opEquals(const T*) const;
    Parameter *param = new Parameter(STCin, type->pointerTo(), Id::p, NULL);
#endif

    parameters->push(param);
    TypeFunction *ftype = new TypeFunction(parameters, Type::tbool, 0, LINKd);
    ftype->mod = MODconst;
    ftype = (TypeFunction *)ftype->semantic(loc, sc);

    FuncDeclaration *fop = new FuncDeclaration(loc, 0, Id::eq, STCundefined, ftype);

    Expression *e = NULL;
    /* Do memberwise compare
     */
    //printf("\tmemberwise compare\n");
    for (size_t i = 0; i < fields.dim; i++)
    {
        Dsymbol *s = fields.tdata()[i];
        VarDeclaration *v = s->isVarDeclaration();
        assert(v && v->storage_class & STCfield);
        if (v->storage_class & STCref)
            assert(0);                  // what should we do with this?
        // this.v == s.v;
        EqualExp *ec = new EqualExp(TOKequal, loc,
            new DotVarExp(loc, new ThisExp(loc), v, 0),
            new DotVarExp(loc, new IdentifierExp(loc, Id::p), v, 0));
        if (e)
            e = new AndAndExp(loc, e, ec);
        else
            e = ec;
    }
    if (!e)
        e = new IntegerExp(loc, 1, Type::tbool);
    fop->fbody = new ReturnStatement(loc, e);

    members->push(fop);
    fop->addMember(sc, this, 1);

    sc = sc->push();
    sc->stc = 0;
    sc->linkage = LINKd;

    fop->semantic(sc);

    sc->pop();

    //printf("-StructDeclaration::buildOpEquals() %s\n", toChars());

    return fop;
}

void CppParticleSimulator::updateSimulation(const Parameters &params, float dt_seconds) {

    // Set forces to 0 and calculate densities
    for (int i = 0; i < positions.size(); ++i) {
        forces[i] = {0, 0, 0};
        float density = 0;

        for (int j = 0; j < positions.size(); ++j) {
            glm::vec3 relativePos = positions[i] - positions[j];
            density += params.get_particle_mass() * Wpoly6(relativePos, params.kernel_size);
        }

        densities[i] = density;
    }

    // Calculate forces
    for (int i = 0; i < positions.size(); ++i) {
        float iPressure = (densities[i] - params.rest_density) * params.k_gas;
        float cs = 0;
        glm::vec3 n = {0, 0, 0};
        float laplacianCs = 0;

        glm::vec3 pressureForce = {0, 0, 0};
        glm::vec3 viscosityForce = {0, 0, 0};

        for (int j = 0; j < positions.size(); ++j) {
            glm::vec3 relativePos = positions[i] - positions[j];

            // Particle j's pressure force on i
            float jPressure = (densities[j] - params.rest_density) * params.k_gas;
            pressureForce = pressureForce - params.get_particle_mass() *
                ((iPressure + jPressure) / (2 * densities[j])) *
                gradWspiky(relativePos, params.kernel_size);

            // Particle j's viscosity force in i
            viscosityForce += params.k_viscosity *
                params.get_particle_mass() * ((velocities[j] - velocities[i]) / densities[j]) *
                laplacianWviscosity(relativePos, params.kernel_size);

            // cs for particle j
            cs += params.get_particle_mass() * (1 / densities[j]) * Wpoly6(relativePos, params.kernel_size);

            // Gradient of cs for particle j
            n += params.get_particle_mass() * (1 / densities[j]) * gradWpoly6(relativePos, params.kernel_size);

            // Laplacian of cs for particle j
            laplacianCs += params.get_particle_mass() * (1 /densities[j]) * laplacianWpoly6(relativePos, params.kernel_size);
        }

        glm::vec3 tensionForce;

        if (glm::length(n) < params.k_threshold) {
            tensionForce = {0, 0, 0};
        } else {
            tensionForce = params.sigma * (- laplacianCs / glm::length(n)) * n;
        }

        //glm::vec3 n = {0, 0, 0};
        glm::vec3 boundaryForce = {0, 0, 0};
        boundaryForce = calculateBoundaryForceGlass(params, i);

        // Add external forces on i
        forces[i] = pressureForce + viscosityForce + tensionForce + params.gravity + boundaryForce;

        // Euler time step
        velocities[i] += (forces[i] / densities[i]) * dt_seconds;
        positions[i] += velocities[i] * dt_seconds;

    }

    checkBoundariesGlass(params);

    glBindBuffer (GL_ARRAY_BUFFER, vbo_pos);
    glBufferData (GL_ARRAY_BUFFER, positions.size() * 3 * sizeof (float), positions.data(), GL_STATIC_DRAW);
}

int PLL::Events (realtype t, N_Vector X, realtype * event, void * data) {
#ifndef WITHPHIERR
  realtype x, y, r, w;
  Parameters * parameters = (Parameters *) data;
  
  x = Ith (X, 0);
  y = Ith (X, 1);
  r = Ith (X, 2);
  w = Ith (X, 3);
  
  fREF = parameters->At(0);
  T = 1.0/fREF;
  R1 = parameters->At(1);
  omega0 = 2*pi*parameters->At(2);
  Vdd = parameters->At(3);
  rho0 = parameters->At(4);
  rhoap = parameters->At(5);
  k0 = parameters->At(6);
  Krho = parameters->At(7);
  Kap = parameters->At(8);
  alpha = parameters->At(9);
  KVCOa = parameters->At(10);
  KVCOb = parameters->At(11);
  KVCOc = parameters->At(12);
  tuning_coeff = parameters->At(13);
  
  // rising edge of the clock
  event[0] = t - n*T;
  // extrema of x
  realtype gamma = sqrt(x*x+y*y);
  realtype vtune = w/tuning_coeff;
  event[1] = ((rho0+Krho*vtune)/gamma - 1)*k0*x - ((1-alpha)*Kap*(gamma-rhoap) + 1 + alpha*vtune*(KVCOa + KVCOb*vtune + KVCOc*vtune*vtune))*omega0*y;
  // reset
  event[2] = t - (treset+tau_d+dt);
#endif
  return CV_SUCCESS;
}

void	Room::setParameters(Parameters &params)
{
  _parameter->setDifficulty(params.getDifficulty());
}

// This is the workhorse of the code.
// It takes in the input parameters and an output object, and runs everything from there.
// This modularization is set up so that anybody can call the evolution routines.
//
// Input parameters:
// inifile: the input parameters
// params: class containing cosmological parameters (somewhat degenerate with inifile, but these values won't be read from inifile)
// output: the outputting class
// postprocess: whether or not to perform postprocessing
//
// Return values:
// 0: success
// -1: error in initialization
// 1: Integration error
// 2: NAN error
// 3: Did not get to a = 1 in allotted time
// 4: Model reports invalid state
// 5: Error in postprocessing integration of distance measures
int doEvolution(IniReader& inifile, Parameters& params, Output& output, vector<vector<double> > &SN1adata, bool postprocess) {

    //****************//
    // Initialization //
    //****************//

    // Set up the integrator
    Integrator myIntegrator;

    // Set up the model class
    Model *myModel;
    std::string parsestring = inifile.getiniString("model", "LambdaCDM", "Cosmology");
    if (parsestring == "Quintessence")
        myModel = new Quintessence();
    else if (parsestring == "LinearW")
        myModel = new LinearW();
    else if (parsestring == "Kessence")
        myModel = new Kessence();
    else if (parsestring == "KGB")
        myModel = new KGB();
    else if (parsestring == "Fr")
        myModel = new Fr();
    else
        myModel = new LambdaCDM();    // LambdaCDM is the default

    // Load the model and parameters into a class to pass into the integration routine
    IntParams myIntParams(params, *myModel);

    // Set up the consistency check class
    parsestring = inifile.getiniString("consistencyclass", "None", "Function");
    Consistency *myChecker;
    if (parsestring == "SimpleCheck")
        myChecker = new SimpleCheck();
    else
        myChecker = new Consistency();  // Default option, which has no checking

    // Create vectors for redshift and hubble
    vector<double> redshift;
    vector<double> hubble;

    int result = 0; // For information coming back from functions


    //********************//
    // Initial Conditions //
    //********************//

    // The initial value of a is extracted from the starting redshift
    // The initial values of \phi and \dot{\phi} are read from the input parameters
    double phi0 = inifile.getiniDouble("phi0", 0.0, "Cosmology");
    double phidot0 = inifile.getiniDouble("phidot0", 0.0, "Cosmology");
    double data[4] = { 1.0 / (1.0 + params.z0()), phi0, phidot0, 0.0 };
    // The data array stores a, \phi, \dot{phi} and H through the evolution
    // H is calculated in the initialization of the model

    // Start and end times
    double starttime = inifile.getiniDouble("starttime", 0.0, "Function");
    double endtime = starttime + inifile.getiniDouble("maxtime", 10.0, "Function");

    // Get the output class to write out information on the run
    output.printinfo(data, params);

    // Allow the model to initialize itself
    result = myModel->init(data, starttime, params, inifile, output);
    if (result != 0) {
        delete myChecker;
        delete myModel;
        return -1;
    }

    // Check that the model has an internally consistent state with the initial data (it should have aborted already if so, but be safe)
    if (myModel->isvalidconfig(data) == false) {
        delete myChecker;
        delete myModel;
        return 4;
    }

    // Write the model name to the output log
    output.printvalue("Model", myModel->classname());
    output.printlog(""); // Whitespace for prettiness


//.........这里部分代码省略.........

FuncDeclaration *StructDeclaration::buildCpCtor(Scope *sc)
{
    //printf("StructDeclaration::buildCpCtor() %s\n", toChars());
    FuncDeclaration *fcp = NULL;

    /* Copy constructor is only necessary if there is a postblit function,
     * otherwise the code generator will just do a bit copy.
     */
    if (postblit)
    {
        //printf("generating cpctor\n");

        StorageClass stc = postblit->storage_class &
                            (STCdisable | STCsafe | STCtrusted | STCsystem | STCpure | STCnothrow);
        if (stc & (STCsafe | STCtrusted))
            stc = stc & ~STCsafe | STCtrusted;

        Parameters *fparams = new Parameters;
        fparams->push(new Parameter(STCref, type->constOf(), Id::p, NULL));
        Type *ftype = new TypeFunction(fparams, Type::tvoid, FALSE, LINKd, stc);
        ftype->mod = MODconst;

        fcp = new FuncDeclaration(loc, 0, Id::cpctor, stc, ftype);

        if (!(fcp->storage_class & STCdisable))
        {
            // Build *this = p;
            Expression *e = new ThisExp(0);
#if !STRUCTTHISREF
            e = new PtrExp(0, e);
#endif
            AssignExp *ea = new AssignExp(0,
                new PtrExp(0, new CastExp(0, new AddrExp(0, e), type->mutableOf()->pointerTo())),
                new PtrExp(0, new CastExp(0, new AddrExp(0, new IdentifierExp(0, Id::p)), type->mutableOf()->pointerTo()))
            );
            ea->op = TOKblit;
            Statement *s = new ExpStatement(0, ea);

            // Build postBlit();
            e = new ThisExp(0);
#if !STRUCTTHISREF
            e = new PtrExp(0, e);
#endif
            e = new PtrExp(0, new CastExp(0, new AddrExp(0, e), type->mutableOf()->pointerTo()));
            e = new DotVarExp(0, e, postblit, 0);
            e = new CallExp(0, e);

            s = new CompoundStatement(0, s, new ExpStatement(0, e));
            fcp->fbody = s;
        }
        else
            fcp->fbody = new ExpStatement(0, (Expression *)NULL);

        members->push(fcp);

        sc = sc->push();
        sc->stc = 0;
        sc->linkage = LINKd;

        fcp->semantic(sc);

        sc->pop();
    }

    return fcp;
}

void createHamiltonian(Geometry const &geometry, DynVars const &dynVars,
                       MyMatrix<std::complex<double> >& matrix,Parameters const &ether,Aux &aux,int type)
// Modified by IS Nov-08-04

{
    int	col, volume, i = 0, p,dir;
    tpem_t hopping,hopping2,bandHop;
    double	 tmp,tmp2;
    int	j,iTmp;
    volume = ether.linSize;
    tpem_t S_ij;
    //static vector<double> phonon_q1(volume);
    //static vector<double> phonon_q2(volume);
    //static vector<double> phonon_q3(volume);


    Phonons<Parameters,Geometry> phonons(ether,geometry);
    for (p = 0; p < matrix.getRank(); p++)
        for (col = 0; col < matrix.getRank(); col++)
            matrix(p,col)=0;

    for (p = 0; p < volume; p++) {
        double phonon_q1=phonons.calcPhonon(p,dynVars,0);
        double phonon_q2=phonons.calcPhonon(p,dynVars,1);
        double phonon_q3=phonons.calcPhonon(p,dynVars,2);
        matrix(p,p) = ether.phononEjt[0]*phonon_q1+ether.phononEjt[2]*phonon_q3+ether.potential[p];
        matrix(p+volume,p+volume) = -ether.phononEjt[2]*phonon_q3+ether.phononEjt[0]*phonon_q1+ether.potential[p];
        matrix(p,p+volume) = (ether.phononEjt[1]*phonon_q2);
        matrix(p+volume,p) = conj(matrix(p,p+volume));

        for (j = 0; j < geometry.z(p); j++) {	/* hopping elements */
            iTmp=geometry.borderId(p,j);;
            if (iTmp>=0) hopping = ether.hoppings[iTmp];
            else hopping= -1.0;
            col = geometry.neighbor(p,j);
            tmp=cos(0.5*dynVars.theta[p])*cos(0.5*dynVars.theta[col]);
            tmp2=sin(0.5*dynVars.theta[p])*sin(0.5*dynVars.theta[col]);
            S_ij=tpem_t(tmp+tmp2*cos(dynVars.phi[p]-dynVars.phi[col]),
                        -tmp2*sin(dynVars.phi[p]-dynVars.phi[col]));

            if (p>col) hopping2 = conj(hopping);
            else hopping2=hopping;
            dir = int(j/2);

            bandHop=ether.bandHoppings[0+0*2+dir*4];
            hopping=hopping2 * bandHop ;
            matrix(p,col) = hopping * S_ij;
            matrix(col,p) = conj(hopping * S_ij);

            bandHop=ether.bandHoppings[0+1*2+dir*4];
            hopping= hopping2 * bandHop;
            matrix(p, col+volume)=hopping * S_ij;
            matrix(col+volume,p)=conj(matrix(p,col+volume));

            //
            bandHop=ether.bandHoppings[1+0*2+dir*4];
            hopping= hopping2 * bandHop;
            matrix(p+volume,col) =  hopping * S_ij;
            matrix(col,p+volume) =  conj(matrix(p+volume,col));

            bandHop=ether.bandHoppings[1+1*2+dir*4];
            hopping=hopping2 * bandHop;
            matrix(p+volume,col+volume) =  hopping * S_ij;
            matrix(col+volume,p+volume) = conj(matrix(p+volume,col+volume));
        }

        if (ether.isSet("tprime")) {
            for (j = 0; j < geometry.z(p,2); j++) {	/* hopping elements */
                iTmp=geometry.borderId(p,j,2);
                hopping= -1.0;
                if (iTmp>=0) hopping = ether.hoppings[iTmp];

                col = geometry.neighbor(p,j,2);
                tmp=cos(0.5*dynVars.theta[p])*cos(0.5*dynVars.theta[col]);
                tmp2=sin(0.5*dynVars.theta[p])*sin(0.5*dynVars.theta[col]);
                S_ij=tpem_t(tmp+tmp2*cos(dynVars.phi[p]-dynVars.phi[col]),
                            -tmp2*sin(dynVars.phi[p]-dynVars.phi[col]));

                if (p>col) hopping2 = conj(hopping);
                else hopping2=hopping;
                dir = int(j/2);

                bandHop=ether.bandHoppings[0+0*2+dir*4]*ether.tprime;
                hopping=hopping2 * bandHop ;
                matrix(p,col) = hopping * S_ij;
                matrix(col,p) = conj(hopping * S_ij);

                bandHop=ether.bandHoppings[0+1*2+dir*4]*ether.tprime;
                hopping= hopping2 * bandHop;
                matrix(p, col+volume)=hopping * S_ij;
                matrix(col+volume,p)=conj(matrix(p,col+volume));

                //
                bandHop=ether.bandHoppings[1+0*2+dir*4]*ether.tprime;
                hopping= hopping2 * bandHop;
                matrix(p+volume,col) =  hopping * S_ij;
                matrix(col,p+volume) =  conj(matrix(p+volume,col));

                bandHop=ether.bandHoppings[1+1*2+dir*4]*ether.tprime;
                hopping=hopping2 * bandHop;
//.........这里部分代码省略.........