scrap.txt

    static void Presence(District *pDistrict, int NTS);
    static void WindowsOpening(Building *pBuilding, Climate *pClimate, unsigned int step, unsigned int step2, double &Nvent);
    static void VentilationNocturne(Building *pBuilding,Climate *pClimate, int step);


///Presence Model developped by Jessen Page in Mathlab and translate in C++ by Christophe Giller
// as input it takes prensence profile for 1 week and long  absence
// The model is explain in more detail in the Jessen Page's Thesis no 3900 (2007).

void Model::Presence(District *pDistrict, int NTS){

    vector<double> pd; //presence distribution
    //load("presence.txt", pd); // commented here JK - 20.03.09 (due to removal of fct.cpp)

    vector<double> T11, T01, T10, vLongAbs,vLongAbs0;
    vector<int> nbPers;
    int nbPersMax;
    double mu, beta, pLongAbsence, x;   //beta is the mobility factor (mu) adjusted
    //load("LongAbs.txt", vLongAbs0); // commented here JK - 20.03.09 (due to removal of fct.cpp)
    for (unsigned int b=0;b<pDistrict->getnBuildings();b++){
    for (unsigned int z=0;z<pDistrict->getBuilding(b)->getnZones();z++){
        mu = pDistrict->getBuilding(b)->getZone(z)->getMu();                   //parametre de mobilité
        pLongAbsence = pDistrict->getBuilding(b)->getZone(z)->getDayOff()/52 / (7*24*4);

        nbPersMax = pDistrict->getBuilding(b)->getZone(z)->getNbPersMax();

        nbPers.assign(NTS, 0);
        T01.assign (pd.size(), 0.0);
        T11.assign (pd.size(), 0.0);
        T10.assign (pd.size(), 0.0);



        ///Determination de T01 et T11
        beta = mu; //ajustement du mu
        for (int i = 0;i < pd.size();i++){
            if(pd[i+1]==0){
                T01[i] = 0;
                T11[i] = 0;
            }
            else if(pd[i+1]==1){
                T01[i] = 1;
                T11[i] = 1;
            }
            else{
                if (pd[i]==1){
                    T11[i] = pd[i+1];
                    T01[i] = 0;
                }
                else if (pd[i]==0){
                    T01[i]=pd[i+1];
                    T11 [i]=0;
                }
                else if (pd[i] == pd[i+1]){
                    if (pd[i] + pd[i+1]>1){
                        if (mu>1/(2*pd[i]-1))
                            beta=1/(2*pd[i]-1);
                        else beta=mu;
                    }


                    else if (pd[i]+pd[i+1] < 1){
                        if (mu > 1/(1-2*pd[i]))
                            beta = 1/(1-2*pd[i]);
                        else beta = mu;
                    }
                    else beta = mu;
                    T01[i] = 2*beta*pd[i]/(beta+1);
                    T11[i] = 1-(1-pd[i])*T01[i]/pd[i];
                }


                else if (pd[i]<pd[i+1]){
                    if (mu<(pd[i+1]-pd[i])/(2-(pd/*chang in pd*/[i+1]+pd/*iciToo*/[i])))
                        beta = (pd[i+1]-pd[i])/(2-(pd/*iciToo*/[i+1]+pd[i]));
                    else if ((pd[i]+pd[i+1] > 1) && (mu>(pd[i]-pd[i+1]+1)/(pd[i+1]+pd[i]-1)))
                            beta=(pd[i]-pd[i+1]+1)/(pd[i+1]+pd[i+1]-1);
                    else if ((pd[i]+pd[i+1]<1) && (mu > (1-pd[i]+pd[i+1])/(1-pd[i]-pd[i+1])))
                            beta=(1-pd[i]+pd[i+1])/(1-pd[i]-pd[i+1]);
                    else beta = mu;

                T01[i] = pd[i+1]+pd[i]*(beta-1)/(beta+1);
                T11[i] = 1/pd[i]*(pd[i+1]-(1-pd[i])*T01[i]);
                }


                else{
                    if (mu<(pd[i]-pd[i+1])/(pd[i+1]+pd[i]))
                        beta =( pd[i]-pd[i+1])/(pd[i+1]+pd[i]);
                    else
                        if ((pd[i]+pd[i+1]>1) && (mu >(pd[i]-pd[i+1]+1)/(pd[i+1]+pd/*change p en pd*/[i]-1)))
                            beta = (pd[i]-pd[i+1]+1)/(pd[i+1]+pd[i]-1);
                        else if ((pd[i]+pd[i+1] < 1) && (mu > (1-pd[i]+pd[i+1])/(1-pd[i]-pd[i+1])))
                            beta=(1-pd[i]+pd[i+1])/(1-pd/*change pd en p*/[i]-pd[i+1]);
                        else beta=mu;

                    T01[i] = pd[i+1]+pd[i]*(beta-1)/(beta+1);
                    T11[i] = 1/pd[i]*(pd[i+1]-(1-pd[i])*T01[i]);
                }
            }
        }

    double dBlockAbs /*duré block absence*/ , absDurationBin = 12.0, absDurationMin =12.0, tEndBlockAbs /*fin du absence block */;

    vLongAbs = vLongAbs0;
    for (int i = 1; i < vLongAbs.size(); i++){
            vLongAbs[i] += vLongAbs[i-1];
            if (vLongAbs[i] > 1.0)
                vLongAbs[i] = 1.0;
    }

    int occ0, occ, blockAbs0, blockAbs, count;

    for (int j = 0;j<nbPersMax;j++){
        occ0=0; //% initial value of occupancy
        blockAbs0=0; //% initial value of block absence: 0 (normal occupancy)
        for (int i = 0;i< NTS-1 ; ){
            for (int k = 0; k < pd.size();k++,i++){
                    if (i >= NTS-1)
                        break;
                    if (blockAbs0==0){ //% currently not in a block absence
                        x=randomUniform(0.0,1.0);
                        if (x < pLongAbsence){
                            blockAbs=1; // start a block absence
                            // determine the duration of block absence
                            x = randomUniform(0.0,1.0);
                            count = 0;
                            for(int l = 0;l < vLongAbs.size();l++){
                                if(vLongAbs[l]<x)
                                    count+=1;
                            }
                            x=randomUniform(0.0,1.0);
                            dBlockAbs=(count+x)*absDurationBin+absDurationMin;      //% hours, minimum absDurationMin
                            tEndBlockAbs = i + dBlockAbs*4;
                        }
                        else blockAbs = 0;

                    }
                    else{ // currently in a block absence, check possible end
                        if (i > tEndBlockAbs) blockAbs=0;
                        else blockAbs=1;
                        }

                    if (blockAbs==1) //% occupant just left on long absence
                        occ=0;
                    else if ((blockAbs0==1)&&(blockAbs==0)) // occupant just came back from long absence
                        occ=1;
                    else{
                        x=randomUniform(0.0,1.0);
                        if (occ0==0){ //% room currently not occupied
                            if (T01[k]>x) occ=1;
                            else occ=0;
                        }
                        else if (occ0==1){// % room currently occupied
                            T10[k]=1-T11[k];
                            if (T10[k]>x) occ=0;
                            else occ=1;
                        }
                    }
                    occ0 = occ;
                    blockAbs0 = blockAbs;

                    nbPers[i+1] += occ;

                }
            }

        }

    pDistrict->getBuilding(b)->getZone(z)->setPresence(nbPers);

    }
    }

}


/// this model was devellopped by Christophe Giller but is not finished, there is still some problem with modeling
/// the mass flow rate thourgh open windows.
void Model::WindowsOpening(Building *pBuilding, Climate *pClimate, unsigned int step, unsigned int step2, double &Nvent){

    double To =  pClimate->getToutCelsius(dt*step), Ti;
    double patm = pClimate->getPatm(step);
    const double C = 0.6; //constante for the mdot in Jessen thesis p 97
    const double g = 9.81;
    double rho; //outside air m V
    double Vdot;
    bool windowOpen = true;
    double winArea, h = 1, W = 1, V; //height of the windows, Volum of the room

    double a, b,c, alea, pactwin;



        for (int i=0;i<pBuilding->getnZones();i++) {
          //  if (step <3840 && step<3840+7*24) presence = 0;
            //else presence = pBuilding->getZone(i)->getPresence(step*12 - 3840 + step2 / 3);


//
//        if (pBuilding->getZone(i)->getPresence(step*4+ step2/3 -Etude*24*4) != 0){
        if (step2==0) Ti = pBuilding->getZone(i)->getTa(step);
        else Ti = pBuilding->getZone(i)->getTaExpl(step*12+step2-1);

//switch (pBuilding->getmodeltype()) {
//  case 1:
//
/////                 1. Simple deterministic model with Ti and To if (Ti > 22 && To < 22 && step >3840) windowOpen = true;
//    case 2 :
/////                 2. Deterministic model with Ti
//                    if (Ti > 25) windowOpen = true;
//    case 3 :
/////                 3. Always closed
//                    windowOpen = false;
//    case 4 :
/////                 4. Always open
//                    windowOpen = true;
//    case 5 :
/////                 5. Elaborated deterministic model with Ti and To
//                   if (Ti > 25 && To < Ti) windowOpen = true;
//    case 6 :
/////                 6. Logistic regression on Ti (Nicol-Humphreys)
//                 a=-10.69; b=0.379;
//                 pactwin = exp(a+b*(Ti))/(1+exp(a+b*(Ti)));
//                 alea = (rand() % 100)/100.0;
//                 if (alea < pactwin) windowOpen = true;
//
//    case 7 :
/////                 7. Logistic regression on To (Nicol-Humphreys)
//                 a=-2.31; b=0.104;
//                 pactwin = exp(a+b*To)/(1+exp(a+b*To));
//                 alea = (rand() % 100)/100.0;
//                 if (alea < pactwin) windowOpen = true;
//    case 8 :
/////                 8. Multiple logistic regression on Ti and To (Rijal et al)
//            a=-6.4; b=0.171; c=0.166;
//            pactwin = exp(a+b*Ti+c*To)/(1+exp(a+b*Ti+c*To));
//            alea = (rand() % 100)/100.0;
//            if (alea < pactwin) windowOpen = true;
//
//
//}
//windowOpen = true;
            if (windowOpen){

 //               winArea = ScalarProduct(pBuilding->getZone(i)->getSurfaceWindow() , pBuilding->getZone(i)->getOpenableRatio());
//                cout<<"winArea "<<winArea<<endl;
//                h = pBuilding->getZone(i)->getHeightWindow(3); //height of the windows
                V = 84;//pBuilding->getZone(i)->getVi(); //Volum of the room

                //Vdot = C/3*winArea*sqrt(abs(Ti-To)*g*h/(((Ti+To)/2.0)+273.15));

                ///modele en tenant compte du vent
//                double windspeed = pClimate->getWindSpeed(step*3600);
//                double anlgeWindDWall = pClimate->getWindDirection(step*3600) - pBuilding->getZone(i)->getAzimuthWall()[3];
//                if (anlgeWindDWall > 90 || anlgeWindDWall < -90) anlgeWindDWall = 90;
//                windspeed = windspeed*cos(anlgeWindDWall*PI/180);
//                double mmair = 28.8/1000;
//                double rhoi = 1.2929 *(273.15/(double)(Ti+273.15))* (patm/101325.0);
                rho = 1.2929 * (patm/101325.0)*(273.15/(double)(To+273.15));
//                cout<< " rhoi "<<rhoi<< " Ti "<<Ti;
//                double pi = rhoi*(Ti+273.15)*8.31/mmair;
//                pi= 101325.0*exp(-(28.97/1000)*9.81*588/*altitude*//(8,3145*(Ti+273.15)));
//                //patm = rho*(To+273.15)*8.31/mmair;
//                double dP = patm - pi + rho*windspeed*windspeed/2;
//                cout<<" dP "<<dP<<endl;
//                Vdot =  C*winArea*sqrt(2*abs(dP)/rho);//*(dP/abs(dP));

//                cout<<"ws "<<windspeed<<"rho "<<rho<<endl;
//                cout << "Patm "<<patm<< " Pi "<<pi<<endl;
                //Vdot=0.038*(winArea/h)*h*sqrt(h)*(sqrt(abs(To-Ti))); //avec la formule a Jessen de matlab
                //with the book from eif W = 1, H =1;
//                double rhoprime = rhoi/(pow(1+(pow(rhoi/rho,0.333)), 3));
//                Vdot = 0.65*0.33*sqrt(8*9.81*1*(rho-rhoi)*rhoprime);

    ///with Roulet
                double mdot = 0.33*rho*h*W*0.6*sqrt((9.81*h*(Ti-To))/(To+273.15));

                Nvent = 3600*Vdot/V; //échange d'air par heure

                //cout<<"mass flow rate "<<mdot<<endl;
                cout << "Window open, air exchange in hour n = " << Nvent << endl;
                cout << "step 1 2 \t"<<step<< " "<<step2<<" Zone "<<i<<endl;

            }
            else Nvent = 0.0;

            pBuilding->getZone(i)->setNvent(Nvent);

        }
       // else Nvent = 0.0;
        //.}
}


void Model::VentilationNocturne(Building *pBuilding,Climate *pClimate, int step){
    double To =  pClimate->getToutCelsius(dt*step), Ti = pBuilding->getZone(0)->getTa(step);;

    if (step % 24 > 22 || step % 24 < 6 && To < Ti && Ti > 16){

        const double C = 0.6; //constante for the mdot in Jessen thesis p 97
        const double g = 9.81;
        double winArea, h = 1, V = pBuilding->getZone(0)->getVi();;

        // JK - commented to remove the ScalarProduct in util.cpp 20.03.2009
        // ScalarProduct is implemented in GENPoint
        //winArea = ScalarProduct(pBuilding->vZone[0]->getSurfaceWindow() , pBuilding->vZone[0]->getOpenableRatio());

        double Vdot = C/3*winArea*sqrt(abs(Ti-To)*g*h/(((Ti+To)/2.0)+273.15));
        double Nvent = 3600*Vdot/V;

        pBuilding->getZone(0)->setNvent(Nvent);
    }
}