Tag: Code

Non classé

Calculation of the incident solar flux on a cylinder

nabladmin

And a routine for calculating the incident solar flux on an individual, represented by a vertical cylinder (method described in this article).

###################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
#  AREP, 16 avenue d'Ivry,75013 Paris, FRANCE
##################################################################
#Calcul du flux recu par un cylindre
#	calcul integral pour le flux direct
#	modelisation anisotrope pour le flux diffus,
#   d'apres la methode de (Perez et. al 1990)
##################################################################
# Last modification : 01/08/2017                                 #
##################################################################
# Copyright (C) 2018 El Mehdi Hamdani et Edouard Walther         #  
# This program is free software; you can redistribute it and/or  #
# modify it under the terms of the GNU General Public License    #
# as published by the Free Software Foundation; either version   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : elmehdi[dot]hamdani[at]arep[dot]fr 
#           edouard[dot]walther[at]arep[dot]fr
##################################################################

import math
global pi
pi = 3.141592653589793

# Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence
def flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence):

	# Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m²
	I0 = 1368
	# angle solaire zenithal
	ksi_z = 90 - math.degrees(hauteur_solaire)
	# composante decrivant la transparence du ciel
	epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3))
	# Masse atmospherique (Kasten and Young 1989)
	a = max(0, math.cos(angle_incidence))
	b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z)))
	ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364))
	

	# Delta represente l'eclairement du ciel
	delta = ma * diffus_horizontal / I0
	
	# coefficients donnes par PEREZ et al (1990)
	if epsilon <= 1.065:
		f11 = 0.013
		f12 = 0.764
		f13 = -0.1
		f21 = -0.058
		f22 = 0.127
		f23 = -0.023
	elif epsilon <= 1.23:
		f11 = 0.095
		f12 = 0.920
		f13 = -0.152
		f21 = 0
		f22 = 0.051
		f23 = -0.02
	elif epsilon <= 1.5:
		f11 = 0.464
		f12 = 0.421
		f13 = -0.28
		f21 = 0.064
		f22 = -0.051
		f23 = -0.002
	elif epsilon <= 1.95:
		f11 = 0.759
		f12 = -0.009
		f13 = -0.373
		f21 = 0.201
		f22 = -0.382
		f23 = 0.01
	elif epsilon <= 2.8:
		f11 = 0.976
		f12 = -0.4
		f13 = -0.436
		f21 = 0.271
		f22 = -0.638
		f23 = 0.051
	elif epsilon <= 4.5:
		f11 = 1.176
		f12 = -1.254
		f13 = -0.462
		f21 = 0.295
		f22 = -0.975
		f23 = 0.129
	elif epsilon <= 6.2:
		f11 = 1.106
		f12 = -1.563
		f13 = -0.398
		f21 = 0.301
		f22 = -1.442
		f23 = 0.212
	else:
		f11 = 0.934
		f12 = -1.501
		f13 = -0.271
		f21 = 0.42
		f22 = -2.917
		f23 = 0.249
	
	# coefficient ponderant les densites du flux circumsolaire
	K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180)
	
	# coefficient ponderant les densites du flux de l'horizon
	K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180
	
	# flux diffus provenant du rayonnement circumsolaire W/m²
	diffus_circumsolaire = diffus_horizontal * K1 * a/b
	
	#retourne la valeur du flux diffus circumsolaire
	return diffus_circumsolaire

def flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence):
# Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence

	# Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m²
	I0 = 1368
	# angle solaire zenithal
	ksi_z = 90 - math.degrees(hauteur_solaire)
	# composante decrivant la transparence du ciel
	epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3))
	
	# Masse atmospherique (Kasten and Young 1989)
	a = max(0, math.cos(angle_incidence))
	b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z)))
	ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364))
	
	# Delta represente l'eclairement du ciel
	delta = ma * diffus_horizontal / I0
	
	# coefficients donnes par PEREZ et al (1990)
	if epsilon <= 1.065:
		f11 = 0.013
		f12 = 0.764
		f13 = -0.1
		f21 = -0.058
		f22 = 0.127
		f23 = -0.023
		
	elif epsilon <= 1.23:
		f11 = 0.095
		f12 = 0.920
		f13 = -0.152
		f21 = 0
		f22 = 0.051
		f23 = -0.02
		
	elif epsilon <= 1.5:
		f11 = 0.464
		f12 = 0.421
		f13 = -0.28
		f21 = 0.064
		f22 = -0.051
		f23 = -0.002
		
	elif epsilon <= 1.95:
		f11 = 0.759
		f12 = -0.009
		f13 = -0.373
		f21 = 0.201
		f22 = -0.382
		f23 = 0.01
		
	elif epsilon <= 2.8:
		f11 = 0.976
		f12 = -0.4
		f13 = -0.436
		f21 = 0.271
		f22 = -0.638
		f23 = 0.051
		
	elif epsilon <= 4.5:
		f11 = 1.176
		f12 = -1.254
		f13 = -0.462
		f21 = 0.295
		f22 = -0.975
		f23 = 0.129
		
	elif epsilon <= 6.2:
		f11 = 1.106
		f12 = -1.563
		f13 = -0.398
		f21 = 0.301
		f22 = -1.442
		f23 = 0.212
		
	else:
		f11 = 0.934
		f12 = -1.501
		f13 = -0.271
		f21 = 0.42
		f22 = -2.917
		f23 = 0.249
		
	# coefficient ponderant les densites du flux circumsolaire
	K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180)
	
	# coefficient ponderant les densites du flux de l'horizon
	K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180
	
	# flux diffus provenant de l'horizon et de la voute celeste W/m²
	diffus_horizon_voute = diffus_horizontal * ((1 - K1) * (0.5 * (1 + math.cos(pi/2))) + K2 * math.sin(pi/2))
	
	# retourne la valeur du flux diffus de l'horizon et de la voute celeste
	return diffus_horizon_voute
	
# Cette fonction permet de retourner la valeur du flux solaire absorbe par un individu, 
# a partir des valeurs des flux solaires diffus et direct horizontaux, de l'albedo du sol, du jour de l'annee, de l'heure, du numero du fuseau, de la latitude et de la longitude du site
def flux_cylindre_anisotrope (diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est): 

	# source : Wikipedia. L'equation du temps permet de calculer le temps solaire vrai.
	equation_temps = 7.678 * math.sin(1.374 + (2 * pi * (jour - 81)/365)) - 9.87 * math.sin(2 * (2 * pi * (jour - 81)/365)) 
	
	# le temps solaire vrai permet de calculer la position du soleil
	temps_solaire_vrai = heure + (longitude_est/15) - fuseau - (equation_temps/60) 
	
	# l'unite de l'angle horaire doit être en radians car dans le langage python l'angle des fonctions trigonometriques est considere comme etant en radians
	angle_horaire = math.radians((temps_solaire_vrai - 12) * 15)
	
	# declinaison en radians
	declinaison = math.asin(0.4 * math.sin(math.radians(0.986 * jour - 80)))
	
	# hauteur solaire en radians
	hauteur_solaire = math.asin(math.sin(math.radians(latitude)) * math.sin(declinaison) + math.cos(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire))
	
	#Afin d'eviter des erreurs dans la suite du calcul nous limitons certaines valeurs de hauteur solaire a une valeur nulle
	if hauteur_solaire < math.radians(2):
		hauteur_solaire = 0
	
	#L'azimut est l'angle compte a partir du sud, positivement vers l'ouest, negativement vers l'est (temps_solaire_vrai < 12), entre le plan vertical du soleil a l'instant donne et le plan meridien local.
	if temps_solaire_vrai < 12:
		azimut = -1 * math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire)) 
	else:
		azimut = math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire))
	
	
	#Cette formule permet de calculer le cosinus de l'angle d'incidence
	#	(angle entre la normale au plan recepteur et le rayon solaire incident)
	# Le math.cos(0) signifie que le plan suit la trajectoire du soleil.
	cosinus_angle_incidence = math.cos(hauteur_solaire)*math.sin(0.5*pi)*math.cos(0)+math.sin(hauteur_solaire)*math.cos(0.5*pi)

	#pour eviter une division par 0 dans la formule ligne 224, lorsque la hauteur solaire est egale a 0 on considere que le flux direct normal est nul	
	if hauteur_solaire == 0:
		direct_normal = 0
	else:
		direct_normal = direct_horizontal/math.sin(hauteur_solaire)

	if diffus_horizontal == 0:
		flux_diffus_circumsolaire = 0
		flux_diffus_voute_horizon = 0	
	else:
		flux_diffus_circumsolaire = flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence))
		# calcul de l'azimut de la face arriere de la paroi verticale en degre
		if math.degrees(azimut) < 0:
			azimut_arr = math.degrees(azimut) + 180
		else:
			azimut_arr = math.degrees(azimut) - 180
		
		azimut_arriere = math.radians(azimut_arr)
		#cosinus de l'angle d'incidence de la face arriere de la paroi verticale
		cosinus_angle_incidence_arriere = math.cos(hauteur_solaire) * math.sin(0.5 * pi) * math.cos(azimut - azimut_arriere) + math.sin(hauteur_solaire) * math.cos(0.5 * pi)
		
		flux_diffus_voute_horizon = flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence)) + flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence_arriere))
	
	if flux_diffus_circumsolaire < 0:
		flux_diffus_circumsolaire = 0

	if flux_diffus_voute_horizon < 0:
		flux_diffus_voute_horizon = 2 * diffus_horizontal * 0.5
	
	if hauteur_solaire == 0:
		flux_direct = 0
	else:
		flux_direct = direct_horizontal * cosinus_angle_incidence / math.sin(hauteur_solaire)
		
	flux_reflechi_direct = albedo * direct_horizontal * cosinus_angle_incidence	* (1 - math.cos(pi/2))/2
	flux_reflechi_diffus = albedo * diffus_horizontal * (1 - math.cos(pi/2))/2
	
	# geometrie du cylindre
	rayon_cylindre = 0.17
	hauteur_cylindre = 1.73
	# calcul des flux recus par le cylindre 
	# 	> direct = anisotrope (calcul integral)
	flux_direct_cylindre = 2 * flux_direct * rayon_cylindre * hauteur_cylindre
	# 	> diffus = anisotrope circumsolaire (calcul integral) + isotrope diffus et horizon
	flux_diffus_cylindre = (2 * flux_diffus_circumsolaire + pi * flux_diffus_voute_horizon) * rayon_cylindre * hauteur_cylindre		
	# 	> reflechi  = anisotrope (calcul integral) + isotrope
	flux_reflechi_cylindre = (2 * flux_reflechi_direct + 2 * pi * flux_reflechi_diffus) * rayon_cylindre * hauteur_cylindre	
		
	#Valeurs en Watts
	return(flux_direct_cylindre, flux_diffus_cylindre, flux_reflechi_cylindre)

def temperature_mrt_fluxsolaire(temperature_mrt_classique, diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est):
	sigma = 5.67 * 0.00000001 #constante de stefan boltzmann
	epsilon = 0.97 #emissivite du corps humain
	#temperature radiante en Kelvin
	temperature_radiante = math.pow(temperature_mrt_classique + 273.15, 4)
	f_eff = 0.75 #facteur de surface de rayonnement effectif pour un individu debout d'après JENDRITZKY and NUBLER
	alpha_clo = 0.8 #coefficient d'absorption des vêtements
	# geometrie cylindre
	rayon_cylindre = 0.17
	hauteur_cylindre = 1.73
	S_cylindre = 2*pi*rayon_cylindre*hauteur_cylindre
	
	#flux direct, diffus, reflechi avec diffus anisotrope en W/m²
	flux_direct, flux_diffus, flux_reflechi = flux_cylindre_anisotrope(diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est) / S_cylindre
	temperature_mrt_solaire = math.pow(temperature_radiante + (alpha_clo * flux_direct + f_eff * alpha_clo * flux_diffus + f_eff * alpha_clo * flux_reflechi)/(sigma * epsilon), 0.25) - 273.15
	
	return temperature_mrt_solaire
Non classé

Scilab / Windows : parallelization of the genetic algorithm optim_GA

nabladmin

Scilab includes a genetic algorithm in the function optim_GA.

To save time when evaluating each individual in a generation, it is possible to parallelize the evaluations: they are independent!

We propose here a modification of the Scilaboptim_ga routine allowing you to evaluate in parallel a number of individuals equal to the number of available threads on the computer (the reduction of computing time is thus proportional to the number of cores of your machine)

This operation is valid under Linux only, the Windows version of Scilab 6 not supporting parallelism: « In this current version of Scilab, parallel_run uses only one core on Windows platforms. » (source).

/////////////////////////////////////////////////////////////////////
//_______  ______    _______  _______ 
//|   _   ||    _ |  |       ||       |
//|  |_|  ||   | ||  |    ___||    _  |
//|       ||   |_||_ |   |___ |   |_| |
//|       ||    __  ||    ___||    ___|
//|   _   ||   |  | ||   |___ |   |    
//|__| |__||___|  |_||_______||___|    
// AREP - 16, av. d'Ivry, 7501/3 Paris, FRANCE
/////////////////////////////////////////////////////////////////////
// Scilab (www.scilab.org) - This file is part of Scilab
// Copyright (C) 2008 - Yann COLLETTE <yann[dot]collette[at]renault[dot]com>
// Copyright (C) 2014 - Michael Baudin <michael[dot]baudin[at]contrib[dot]scilab[dot]org>
//
//	26/6/2017 : added the "parallel_run" capability for the evaluation of individuals (Edouard Walther)
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt                                  //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Last modification : 26/06/2017                                 //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Contact : edouard[dot]walther[at]arep[dot]com                                
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_GA_parallel(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param)

    [nargout, nargin] = argn();

    if ~isdef("param", "local") then
        param = [];
    end

    [codage_func, err]    = get_param(param, "codage_func", coding_ga_identity);
    [init_func, err]      = get_param(param, "init_func", init_ga_default);
    [crossover_func, err] = get_param(param, "crossover_func", crossover_ga_default);
    [mutation_func, err]  = get_param(param, "mutation_func", mutation_ga_default);
    [selection_func, err] = get_param(param, "selection_func", selection_ga_elitist);
    [nb_couples, err]     = get_param(param, "nb_couples", 100);
    [pressure, err]       = get_param(param, "pressure", 0.05);
    [output_func, err] = get_param(param, "output_func", output_ga_default);

    if ~isdef("ga_f", "local") then
        error(sprintf(gettext("%s: ga_f is mandatory"), "optim_ga"));
    else
        if typeof(ga_f) == "list" then
            deff("y = _ga_f(x)", "y = ga_f(1)(x, ga_f(2:$))");
        else
            deff("y = _ga_f(x)", "y = ga_f(x)");
        end
    end

    if ~isdef("pop_size", "local") then
        pop_size = 100;
    end
    if ~isdef("nb_generation", "local") then
        nb_generation = 10;
    end
    if ~isdef("p_mut", "local") then
        p_mut = 0.1;
    end
    if ~isdef("p_cross", "local") then
        p_cross = 0.7;
    end
    if ~isdef("Log", "local") then
        Log = %F;
    end

    // Initialization of the population
    Pop = list();
    Pop = init_func(pop_size, param);

    if (nargout >= 3) then
        pop_init = Pop;
    end

    // Code the individuals
    Pop = codage_func(Pop, "code", param);

	// Getting the objective function for each individual
	//disp("Extraction list...");
	[vec_param_pop]=Pop(1);
	vec_param_total=vec_param_pop;
    for i = 2:length(Pop)
        [vec_param_pop]=Pop(i);
		vec_param_total=cat(2,vec_param_total,vec_param_pop);
    end

	// First launch 
	//disp("Execution for the first population...");
	[FObj_Pop]=parallel_run(vec_param_total,_ga_f);
	FObj_Pop=FObj_Pop';
		
    if (nargout == 4) then
        fobj_pop_init = FObj_Pop;
    end

    FObj_Pop_Max = max(FObj_Pop);
    FObj_Pop_Min = min(FObj_Pop);

    // Normalization of the efficiency
    Efficiency = (1 - pressure) * (FObj_Pop_Max - FObj_Pop) / max([FObj_Pop_Max - FObj_Pop_Min %eps]) + pressure;

	
	//disp("Starting optimisation with GA...");
	
    // The genetic algorithm
    for i = 1:nb_generation
        //
        // Selection
        //
        Indiv1 = list();
        Indiv2 = list();
        Wheel = cumsum(Efficiency);
        for j = 1:nb_couples
            // Selection of the first individual in the couple
            Shoot = grand(1, 1, "unf", 0, Wheel($));
            Index = find(Shoot <= Wheel, 1);
            Indiv1(j)      = Pop(Index);
            FObj_Indiv1(j) = FObj_Pop(Index);
            // Selection of the second individual in the couple
            Shoot = grand(1, 1, "unf", 0, Wheel($));
            Index = 1;
            Index = find(Shoot <= Wheel, 1);
            Indiv2(j)      = Pop(Index);
            FObj_Indiv2(j) = FObj_Pop(Index);
        end
        //
        // Crossover
        //
        for j = 1:nb_couples
            if (p_cross>grand(1, 1, "def")) then
                [x1, x2] = crossover_func(Indiv1(j), Indiv2(j), param);
                Indiv1(j) = x1;
                Indiv2(j) = x2;
                ToCompute_I1(j) = %T;
                ToCompute_I2(j) = %T;
            else
                ToCompute_I1(j) = %F;
                ToCompute_I2(j) = %F;
            end
        end
        //
        // Mutation
        //
        for j = 1:nb_couples
            if (p_mut>grand(1, 1, "def")) then
                x1 = mutation_func(Indiv1(j), param);
                Indiv1(j) = x1;
                ToCompute_I1(j) = %T;
            end
            if (p_mut>grand(1, 1, "def")) then
                x2 = mutation_func(Indiv2(j), param);
                Indiv2(j) = x2;
                ToCompute_I2(j) = %T;
            end
        end
        //
        // Computation of the objective functions
		
		k=0;kk=0; // counters to iterate 
		for j = 1:nb_couples // for all couples in the population
			if ToCompute_I1(j) then// if to be computed
				k=k+1;
				if k==1 then // create the first vector of parameters
					[vec_param_pop1]=Indiv1(j);
				else // concatenate for parallel_run
					[vec_param_indiv1]=Indiv1(j);
					indices_indiv1(k)=j;
					vec_param_pop1=cat(2,vec_param_pop1,vec_param_indiv1);
				end
			end
			if ToCompute_I2(j) then// if to be computed
				kk=kk+1;
				if kk==1 then
					[vec_param_pop2]=Indiv2(j);
				else
					[vec_param_indiv2]=Indiv2(j);
					indices_indiv2(kk)=j;
					vec_param_pop2=cat(2,vec_param_pop2,vec_param_indiv2);
				end
			end
		end
		
		// Parallel_run
		//disp("Parallel launch for Indiv1...");
		[objectifs_Indiv1]=parallel_run(vec_param_pop1,_ga_f);
		objectifs_Indiv1=objectifs_Indiv1';
		//disp("Parallel launch for Indiv2...");
		[objectifs_Indiv2]=parallel_run(vec_param_pop2,_ga_f);
		objectifs_Indiv2=objectifs_Indiv2';
		
		// Updating indexes
		//disp("Updating FObj1 ...");
		for k=1:length(objectifs_Indiv1)
			if indices_indiv1(k)<> 0 then
				FObj_Indiv1(indices_indiv1(k))= objectifs_Indiv1(k);
			end;
		end
		for k=1:length(objectifs_Indiv2)
			if indices_indiv2(k)<> 0 then
				FObj_Indiv2(indices_indiv2(k))= objectifs_Indiv2(k);
			end;
		end

        // Reinit ToCompute lists
        ToCompute_I1 = ToCompute_I1 & %F;
        ToCompute_I2 = ToCompute_I2 & %F;
        // Recombination
        [Pop, FObj_Pop] = selection_func(Pop, Indiv1, Indiv2, FObj_Pop, FObj_Indiv1, FObj_Indiv2, [], [], [], param);
        // Callback for plotting / printing intermediate results or stopping the algorithm
        if (Log) then
            stop = output_func(i, nb_generation, Pop, FObj_Pop, param);
            if (stop) then
                break
            end
        end
    end

    pop_opt  = Pop;
    pop_opt  = codage_func(pop_opt, "decode", param);
    fobj_pop_opt = FObj_Pop;
endfunction

 

Non classé

S.E.T. calculation in Python

nabladmin

Still in the comfort models, the SET* comfort index code is given below. The index is calculated according to the most up-to-date models described in the article "Comfort Modeling" on this website.

###################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
#  AREP, 16 avenue d'Ivry,75013 Paris, FRANCE
##################################################################
#Calcul de la SET en Python
##################################################################
# Last modification : 12/03/2018                                 #
##################################################################
# Copyright (C) 2018 Edouard Walther                             #  
# This program is free software; you can redistribute it and/or  #
# modify it under the terms of the GNU General Public License    #
# as published by the Free Software Foundation; either version   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]fr
##################################################################

import math
import random

def pv_sat(T):
	if T >= 0:
		pv_sat2 = pow(10, (2.7877 + (7.625*T)/(241 + T)))
	else :
		pv_sat2 = pow(10,(2.7877 + (9.756*T)/(272.7 + T)))
	return pv_sat2

def pv_calc(T, RH):
	pv_calc2 = (RH * pv_sat(T))/100
	return pv_calc2

def w(T, RH, p):
	w2 = 0.622 * (pv_calc(T, RH)/(p - pv_calc(T, RH)))
	return w2

def Cp_ah(T, RH, p):
	cpa = 1006
	cpv = 1830
	water = w(T, RH, p)
	Cp_ah2 = (cpa + water * cpv)/(1 + water)
	return Cp_ah2

def v_spe(T, RH, p):
	v_spe2 = (461.24 * (T + 273.15) * (0.622 + w(T, RH, p)))/p
	return v_spe2

def h_radiation(T_rad, T_surf):
	h_radiation2 = 0.72 * 5.67 * 0.00000001 * ((T_rad + T_surf) + 2 * 273.15) * (pow((T_rad + 273.15), 2) + pow((T_surf + 273.15), 2)) * 0.97 #<--- emissivite = 0,97
	return h_radiation2


###############################################################################
############################ FONCTION SET* ####################################
############################################################################### 	

def modele_metabolique_SET(RadTempMtx, WindSpeedMtx, T_ambient, phi_ambient, p_ambient, hauteur, masse, fat, Cst_dilat,Cst_sweat, Cst_constr, T_core_set, T_skin_set, SkinBloodFlow_set,U_muscle_fat_skin, C_shiv):
	T_skin = T_skin_set
	T_core = T_core_set
	dT = 60
	#p_ambient = 101325
	#met = 1.1
	#clo = 1
	#i_m = 0.45
	# metabolisme masculin en W
	age=30
	surface =  0.203*pow(hauteur,0.725)*pow(masse,0.425) #Surface exterieure du sujet [m2] 
	#genre=random.random()
	genre = 0.1 # homme
	if genre<0.5:
		metabolisme_W = 3.45 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.01 * (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 43.4))
	else:
		metabolisme_W = 3.19 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.018* (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 42.1))
	met = metabolisme_W/surface/58.2
	#print surface, met
	SkinBloodFlow = 6.3
	minutes_metab = 60
	minutes = minutes_metab - 1	
	duree=(minutes + 1) * 60
	temps=0
	while temps < duree: 
		#variables fonctionnelles 
		compteur = 0 #variable de boucle
		temps = temps + dT
		i_m = 0.45
		i_m_static = i_m
		clo = 1 #0.155 m2.K/W
		clo_static = clo
		#evolution metabolique dynamique
		#if temps < duree/2:
		#	met = 2.2
		#	v_walk = 4/3.6
		#else:
		#	met = 1.2
		#	v_walk = 0  	
		  	#WindSpeedMtx = 0.2	
		  	#clo dynamique
		#if v_walk < 0.7:
		#	v_walk = 0.0052 * (met * 58.2 - 58)
		#corr_T = math.exp(0.042 - 0.398 * WindSpeedMtx + 0.066 * WindSpeedMtx**2 - 0.378 * v_walk + 0.094 * v_walk**2)	
		#if WindSpeedMtx > 3.5:
		#	corr_T = 0.582
		# ATTENTION : ICI CHANGER "WindSpeedMtx" car sinon on le change a chaque pas de temps !
		#WindSpeedMtx = math.sqrt(WindSpeedMtx**2 + v_walk**2)
		#clo = clo_static * corr_T
		#i_m = i_m_static * (4.9 - 6.5 * corr_T + 2.6 * corr_T**2)
		#Constantes du corps humain
		KCLO = 0.25  #coefficient augmentation surface d'echange
#		masse = vect[1] #masse moyenne sur une population
#		R_muscle_fat_skin = 5.28
		#Cp_body = 0.97*3600 #capacite calorifique du corps humain [J/(kg.K)]
		body_mass = masse
#		fat = 15 # pourcentage masse graisseuse
		fat_mass = fat/100*body_mass
		Cp_body = fat_mass/body_mass*2510 + (body_mass - fat_mass)/body_mass*3650 # Modele HAVENITH
	  	#Constantes de regulation de l'organisme
		SBC = 0.0000000567 #Constante de Stefann-Boltzmann
#		Cst_sweat = 170
#		Cst_dilat = 200
#		Cst_constr = 0.5
		#Valeurs de consignes de la regulation du corps humain
#		T_skin_set = 33.7 #temperature de peau
#		T_core_set = 36.8 #temperature interne
		T_body_set = 0.1*T_skin_set + 0.9*T_core_set #temperature corporelle
#		SkinBloodFlow_set = 6.3 #debit sanguin [L/m2.h]
		#Conversion d'unites
		p_atmosphere = p_ambient/1000 #conversion en kPa
		p_atmosphere = p_atmosphere*0.009869 #conversion en atm
		R_clo = 0.155*clo
		#correction de la veture  ########## 
		if clo < 0.5:
			f_surf_clo = 1 + 0.2*clo
			
		else:
			f_surf_clo = 1 + 0.15*clo
		
		#calcul du nombre de Lewis
		Lewis = 2434 * v_spe(T_ambient, phi_ambient*100, p_ambient)/(Cp_ah(T_ambient, phi_ambient*100, p_ambient) * 1.04 * pow(0.83, (2/3))) * (18/8.32/(T_ambient + 273.15))
		#Calculs initiaux du metabolisme
		RM = met * 58.2
		Metab = met * 58.2
		w_crit = 0.59 * pow(WindSpeedMtx, (-0.08))
		#Calcul des coefficients d'echange convectif
		h_c = 3 * pow(p_atmosphere, 0.53)
		h_c_vent = 8.600001 * pow(WindSpeedMtx * p_atmosphere, 0.53)
		h_c = max(h_c, h_c_vent)
		#Coefficient d'echange radiatif
		h_r = 4.7
		#Coefficient d'echange global
		h_g = h_r + h_c
		#Resistance thermique convective+radiative
		R_air = 1/(f_surf_clo * h_g)
		#Temperature operative
		T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g #Sous forme de matrice		
		#Temperature superficielle de veture
		T_clo = T_op + (T_skin - T_op)/(h_g * (R_air + R_clo))
		T_clo_OLD = T_clo + 0.5	      
	############### Boucle calcul T_clo ###########################################
		while abs(T_clo - T_clo_OLD) > 0.001:
	  		T_clo_OLD = T_clo
	  		h_r = h_radiation(RadTempMtx, T_clo)# Sous forme de matrice
	  		h_g = h_r + h_c
	  		R_air = 1/(f_surf_clo * h_g)
	  		T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g# Sous forme de matrice
	  		T_clo = (R_air * T_skin + R_clo * T_op)/(R_air + R_clo)
	  		compteur = compteur + 1
	  		if compteur > 20:
	  			break
	###############################################################################
	  	#SkinBloodFlow = SkinBloodFlow_set
	  	#Temperature corporelle
	  	#alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
	  	#T_body = alpha * T_skin + (1 - alpha) * T_core	  			
	################ Modele de regulation du corps humain #######################
		#Skin signal
		signal_skin = T_skin - T_skin_set
		if signal_skin > 0:
			warm_skin = signal_skin
			cold_skin = 0
		else :
	  		warm_skin = 0
	  		cold_skin = -signal_skin
		# Core signal
		signal_core = T_core - T_core_set
		if signal_core > 0:
			warm_core = signal_core
			cold_core = 0
		else:
			warm_core = 0
			cold_core = -signal_core
		#Debit sanguin
		SkinBloodFlow = (SkinBloodFlow_set + Cst_dilat * warm_core)/(1 + Cst_constr * cold_skin)
		if SkinBloodFlow > 90 :
			SkinBloodFlow = 90
		if SkinBloodFlow < 0.5 :
			SkinBloodFlow = 0.5
		#Temperature corporelle
		alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
		T_body = alpha * T_skin + (1 - alpha) * T_core
		#Corps/Body
		signal_body = T_body - T_body_set
		if signal_body > 0 :
			warm_body = signal_body
			cold_body = 0
		else :
			warm_body = 0
			cold_body = -signal_body
		#Debit sudation
		qm_sweat = Cst_sweat * warm_body * math.exp(warm_skin/10.7)
		if qm_sweat > 500:
			qm_sweat = 500
		#Chaleur latente maximale echangee
		R_vap_tot = (R_clo + R_air)/(Lewis * i_m)
		E_max = (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))/R_vap_tot
		h_e = (2.2 * h_c)/(1 + 0.928 * R_clo * h_c)/133.322
		#E_max = h_e * (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))
		#Chaleur latente sudation
		E_sweat = 0.68 * qm_sweat
		#Ratios
		pcent_sweat = E_sweat/E_max #Part d'energie echangee due a la sudation adimensionnel
		pcent_wet = 0.06 + 0.94 * pcent_sweat #Part de la surface du corps mouille
		#Chaleur latente echangee par diffusion de l'eau a travers la couche cutanee
		E_diff = pcent_wet * E_max - E_sweat
		#Chaleur latente totale echangee par la peau
		E_skin = E_sweat + E_diff
		#Sudation supercritique
		if pcent_wet > w_crit :
	  		pcent_wet = w_crit
	  		pcent_sweat = w_crit/0.94
	  		E_sweat = pcent_sweat * E_max
	  		E_diff = 0.06 * (1 - pcent_sweat) * E_max
	  		E_skin = E_sweat + E_diff
	  		drip_cond_nope = 1
		#Condensation (la pression de vapeur saturante a la temperature de la peau est inferieure a la pression de vapeur de l'air ambiant
		elif E_max < 0 :
	  		E_diff = 0
	  		E_sweat = 0
	  		pcent_wet = w_crit
	  		pcent_sweat = w_crit
	  		E_skin = E_max
	  		drip_cond_nope = -2
		else:
			drip_cond_nope = 0
		w_skin = pcent_wet
		#Frisson/shivering
		M_shiv = C_shiv * cold_skin * cold_core
		Metab = RM + M_shiv
		w_sweat_global = E_sweat/E_skin #Part d'humidite due a la sudation
		#Flux echange entre l'interieur du corps (noyau) et la peau
		Flx_core_skin = (T_core - T_skin) * (U_muscle_fat_skin + 1.163 * SkinBloodFlow) # 1.163 = 4 200/3 600 * 1000
		#Chaleur sensible echangee entre la peau et l'exterieur
		DRY = (T_skin - T_op)/(R_air + R_clo)
		#Chaleur sensible echangee par la respiration
		T_expir = 32.5 + 0.066 * T_ambient + 1.99 * 0.000001 * phi_ambient * pv_sat(T_ambient)
		C_resp = 0.0014 * Metab * (T_expir - T_ambient) #
		#Chaleur latente echangee par la respiration
		E_resp = 0.000017251 * Metab * (pv_sat(35.5) - phi_ambient * pv_sat(T_ambient))
		
		#Accumulation d'energie par la peau
		SSK = Flx_core_skin - DRY - E_skin
		Accumulation_skin = SSK
		
		#Accumulation d'energie par le corps
		SCR = Metab - Flx_core_skin - E_resp - C_resp
		Accumulation_core = SCR
		
		#Modification des temperatures par l'effet de l'accumulation
		
		#Methode 1
		dT_skin = SSK * surface * dT /(alpha * masse * Cp_body)
		dT_core = SCR * surface * dT / ((1-alpha) * masse * Cp_body)
		
		#Methode 2
		TCSK = Cp_body/3600 * alpha * masse
		TCCR = Cp_body/3600 * (1 - alpha) * masse
		
		DTSK = (SSK * surface)/(TCSK * dT)
		DTCR = (SCR * surface)/(TCCR * dT)
	
		T_skin = T_skin + DTSK
		T_core = T_core + DTCR
		T_body = alpha * T_skin + (1 - alpha) * T_core		
		#Energie totale echangee par la peau
		H_skin = DRY + E_skin
		#Metabolisme net
		RN = Metab
		#Chaleur latente echangee par sudation en etat de confort
		E_comf = 0.42 * (RN - (1*58.2))
		if E_comf < 0:
			E_comf = 0
		E_max = E_max * w_crit
		#Chaleur latente evaporative requise pour la thermoregulation
		E_req = RN - E_resp - C_resp - DRY
		E_sweat_global = E_sweat
		#Coefficient d'echange chaleur sensible
		HD = 1/(R_air + R_clo)
		#Coefficient d'echange evaporatif
		HE = 1/R_vap_tot
		#Pression de vapeur saturante a la temperature cutanee
		PSSK = pv_sat(T_skin)
		#Coefficients d'echange
		h_r_SET = h_r
		if met < 0.85:
			h_c_SET = 3
		else:
			h_c_SET = 5.66 * pow((met - 0.85), 0.39)
			if h_c_SET < 3:
				h_c_SET = 3
		h_g_SET = h_c_SET + h_r_SET
		#CLO metabolique
		RCLOS = 1.52/(met + 0.6944) - 0.1835
		#Resistance thermique de la veture
		RCLS = 0.155 * RCLOS
		#Correction de la surface d'echange due a la veture
		f_surf_clo_SET = 1 + KCLO * RCLOS
		#Facteur d'efficacite de BURTON
		F_clo_SET = 1/(1 + 0.155 * f_surf_clo_SET * h_g_SET * RCLOS)
		#Index de permeabilite de la veture
		i_m_SET = 0.38
		#Indice de permeation de la veture
		i_clo_SET = i_m_SET * h_c_SET/h_g_SET * (1 - F_clo_SET)/(h_c_SET/h_g_SET - F_clo_SET * i_m_SET)
		#Resistance convective + radiative corrigee
		R_air_SET = 1/(f_surf_clo_SET * h_g_SET)
		#Resistances evaporatives
		REAS = 1/(Lewis * f_surf_clo_SET * h_c_SET)
		RECLS = RCLS/(Lewis * i_clo_SET)
		#Resistance totale au transfert de chaleur sensible
		HD_S = 1/(R_air_SET + RCLS)
		HE_S = 1/(REAS + RECLS)
		#Variables de resolution de SET/ET
		Delta = 0.0001
		dx = 100
		X_OLD = T_skin - H_skin/HD_S
		while abs(dx) > 0.001:
			#compteurr=compteurr+1
			ERR1 = H_skin - HD_S * (T_skin - X_OLD) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD))
			ERR2 = H_skin - HD_S * (T_skin - (X_OLD + Delta)) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD + Delta))
			#if ERR2==ERR1:
			#	print("attention sortie brutale de boucle...")
			#	break
			x = X_OLD - Delta * ERR1/(ERR2 - ERR1)
			dx = x - X_OLD
			X_OLD = x
		ET_global = x
	return ET_global
Non classé

The corrected P.E.T.

nabladmin

Here is the modified P.E.T. comfort index, as described in thefollowing article. The modifications are taken into account.

##################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
# AREP - 16, av. d'Ivry, 75013 Paris, FRANCE
##################################################################
# based on: Peter Hoeppe PET fortran code, from:
# "Urban climatic map and standards for wind environment - Feasibility study, Technical Input Report No.1",
# Edouard Walther and Quentin Goestchel
# Most of the changes in the implementaion are explained in the resolution function comments                                  #
##################################################################
# Last modification : 10/04/2018                                 #
##################################################################
# Copyright (C) 2018 Édouard WALTHER                             #  
# This program is free software; you can redistribute it and/or  #
# modify it under the terms of the GNU General Public License    #
# as published by the Free Software Foundation; either version   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]com                                
##################################################################


import os
import numpy as np
import math as math
import scipy.optimize as optimize

## Implementation of the skin and core temperatures set values #######
tc_set=36.6 # 36.8
tsk_set=34 # 33.7
tbody_set=0.1*tsk_set+0.9*tc_set # Calculation of the body temperature through a weighted average


## Skin blood flow calculation function: #######
def vasoC(tcore,tsk):
    #Set value signals to consider every cases:
    sig_skin = tsk_set - tsk
    sig_core = tcore - tc_set
    if sig_core<0:
        # In this case, Tcore<Tc_set: the body needs to keep the heat --> the blood flow is reduced
        sig_core=0.
    if sig_skin<0:
        # In this case, Tsk>Tsk_set: the body needs to loose heat --> the blood flow is increased
        sig_skin=0.
    # 6.3 L/m^2/h is the set value of the blood flow
    qmblood = (6.3 + 75. * sig_core) / (1. + 0.5 * sig_skin)
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmblood>90:
        qmblood=90.
    # Alpha can be used to calculate tbody in ameliorated models
    #alfa = 0.04177 + 0.74518 / (qmblood + 0.585417)
    alfa = 0.1
    return (qmblood,alfa)


## Sweating flow calculation function: #######
def Suda(tbody,tsk):
    sig_body = tbody - tbody_set
    sig_skin = tsk - tsk_set
    if sig_body<0:
        #In this case, Tbody<Tbody_set: the body needs to keep the heat --> The sweat flow is 0
        sig_body=0.
    if sig_skin<0:
        # In this case, Tsk<Tsk_set: the body needs to keep the heat --> the sweat flow is reduced
        sig_skin=0.
    #qmsw = 170 * sig_body * math.exp((sig_skin) / 10.7)  # [g/m2/h] Expression from Gagge Model
    qmsw = 304.94*10**(-3) * sig_body
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmsw > 500:
        qmsw = 500
    return (qmsw)


## Vectorial MEMI balance calculation function for the 3 node model: #######
def Syst(T, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,mode):
    ## Conversion of T vector in an array
    arr = np.ones((3,1))
    arr[0,0]=T[0] #Corresponds to T_core
    arr[1,0]=T[1] #Corresponds to T_skin
    arr[2,0]=T[2] #Corresponds to T_clothes
    T=arr
    enbal_vec = np.zeros((3,1)) #Useful for the vectorial expression of the balance

    ## Area parameters of the body: #######
    Adu = 0.203 * mbody ** 0.425 * ht ** 0.725
    feff=0.725
    if pos == 1 or pos == 3:
        feff = 0.725
    if pos == 2:
        feff = 0.696
    # Calculation of the Burton coefficient, k = 0.31 for Hoeppe:
    fcl = 1 + (0.31 * icl) # Increase of the heat exchange surface  depending on clothing level
    facl = (173.51 * icl - 2.36 - 100.76 * icl * icl + 19.28 * icl ** 3.0) / 100
    Aclo = Adu * facl + Adu * (fcl - 1.0)
    Aeffr = Adu * feff  # Effective radiative area depending on the position of the subject
    # Partial pressure of water in the air depending on relative humidity and air temperature:
    if mode: # actual environment
        vpa = HR / 100.0 * 6.105 * math.exp(17.27 * Ta / (237.7 + Ta )) #[hPa]
    else:# mode==False implies we are calculating the PET
        vpa= 12 # [hPa] vapour pressure of the standard environment

    # Convection coefficient depending on wind velocity and subject position
    hc = 0
    if pos == 1:
        hc = 2.67 + (6.5 *v**0.67)
    if pos == 2:
        hc = 2.26 + (7.42 *v**0.67)
    if pos == 3:
        hc = 8.6 * (v ** 0.513)
    hc = hc * (p / po) ** 0.55

    # Basic metabolism for men and women in [W] #######
    metab_female = 3.19 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.018 * (ht * 100.0 / mbody**(1.0/3.0) - 42.1))
    metab_male = 3.45 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.01 * (ht * 100.0 / mbody**(1.0/3.0) - 43.4))
    # Source term : metabolic activity
    eswpot = (M + metab_male)/Adu
    fec = (M + metab_female)/Adu
    he = 0.0
    # Attribution of internal energy depending on the sex of the subject
    if sex == 1:
        he = eswpot
    elif sex == 2:
        he = fec
    h = he *(1.0 - eta)  # [W/m2]

    # Respiratory energy losses
    # Expired air temperature calculation:
    texp = 0.47 * Ta + 21.0  # [degC]

    # Pulmonary flow rate
    dventpulm = he * 1.44 * 10.0**(-6.0)

    # Sensible heat energy loss:
    eres = cair * (Ta - texp) * dventpulm  # [W/m2]

    # Latent heat energy loss:
    vpexp = 6.11 * 10.0**(7.45 * texp / (235.0 + texp))
    erel = 0.623 * Lvap / p * (vpa-vpexp) * dventpulm  # [W/m2]
    ere = eres + erel  # [W/m2]

    # Clothed fraction of the body approximation
    rcl = icl / 6.45  # Conversion in m2.K/W
    y = 0
    if facl > 1.0:
        facl = 1.0
    if icl >= 2.0:
        y = 1.0
    if icl > 0.6 and icl < 2.0:
        y = (ht - 0.2)/ht
    if icl <= 0.6 and icl > 0.3:
        y = 0.5
    if icl <= 0.3 and icl > 0.0:
        y = 0.1
    # calculation of the closing radius depending on the clothing level (6.28 = 2* pi !)
    r2 = Adu * (fcl - 1.0 + facl) / (6.28 * ht * y)  # External radius
    r1 = facl * Adu /(6.28 * ht * y)  # Internal radius
    di = r2 - r1

    # Calculation of the thermal resistance of the body
    alpha = vasoC(T[0,0],T[1,0])[1]
    tbody = alpha * T[1,0] + (1 - alpha) * T[0,0]
    htcl = (6.28 * ht * y * di)/(rcl * math.log(r2 / r1)*Aclo)  # [W/(m2.K)]

    # Calculation of sweat losses

    qmsw = Suda(tbody,T[1,0])
    # Lvap Latent heat of evaporation : 2400 [J/g] divided by 3600 for [g/m2/h] to [g/m2/s]
    esw = 2400 * qmsw / 3600  # [W/m2]
    # Saturation vapor pressure at temperature Tsk and for 100% HR
    Pvsk = 6.105 * math.exp((17.27 * (T[1,0]+273.15) - 4717.03)/ (237.7 + T[1,0])) # [hPa]

    rscl=0.155*icl
    Lw = 16.7 * 10 ** (-1)  # [K/hPa] Lewis factor
    he_diff = hc * Lw
    fecl=1/(1+0.92*hc*rscl)
    emax = he_diff * fecl * (Pvsk - vpa)
    w = esw / emax  # Dermal wetness
    if w > 1:
        w=1
        delta = esw-emax
        if delta < 0:
            esw=emax
    if esw < 0:
        esw=0
    i_m=0.3
    R_ecl=(1/(fcl*hc) + rscl)/(Lw*i_m)
    #R_ecl=0.79*1e7 #version Hoeppe
    ediff = (1 - w)*(Pvsk - vpa)/R_ecl  # Basal perspiration
    evap = -(ediff + esw)  # [W/m2]

    # Radiation losses
    # For bare skin area:
    rbare = Aeffr*(1.0 - facl) * emsk * sigm * ((Tmrt + 273.15)**(4.0) - (T[1,0] + 273.15)**(4.0))/Adu
    # For dressed area:
    rclo = feff * Aclo * emcl * sigm * ((Tmrt + 273.15)**(4.0) - (T[2,0] + 273.15)**(4.0))/Adu
    rsum = rclo+rbare


    ## Convection losses #######
    cbare = hc * (Ta - T[1,0]) * Adu * (1.0 - facl)/Adu  # [w/m^2]
    cclo = hc * (Ta - T[2,0]) * Aclo/Adu  # [W/m^2]
    csum = cclo+cbare

    ## Balance equations of the 3-nodes model
    enbal_vec[0,0] = h + ere - (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) # Core balance [W/m^2]
    enbal_vec[1,0] = rbare + cbare + evap + (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) - htcl*(T[1,0]-T[2,0])  # Skin balance [W/m^2]
    enbal_vec[2,0] = cclo + rclo + htcl*(T[1,0]-T[2,0]) # Clothes balance [W/m^2]
    enbal_scal = h + ere + rsum + csum +evap

    if mode:
        return [enbal_vec[0,0],enbal_vec[1,0],enbal_vec[2,0]] #List of the balances values
    else:
        return enbal_scal #Scalar balance used in the PET calculation

# Solving the system
def resolution(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl, Tx):
    Tn = optimize.fsolve(Syst,Tx ,args=(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,True))
    return (Tn, 1)


# PET calculation with dichotomy method 
def PET (age, sex, ht, mbody, pos, M, icl, Tstable,a,b,eps):
    # Definition of a function with the input variables of the PET reference situation
    def f(Tx):
        return Syst(Tstable, Tx, Tx, 50, 0.1, age, sex, ht, mbody, pos, M, 0.9,False)
    Ti = a # Start index of the browsing interval
    Tf = b # End index of the browsing interval
    pet = 0
    while Tf-Ti>eps: # Dichotomy loop
        if f(Ti)*f(pet)<0:
            Tf = pet
        else:
            Ti = pet
        pet = (Ti + Tf) / 2.0
    return pet

# Input data
# Environment constants
po = 1013.25 #[hPa]
rob = 1.06 # Blood density kg/L
cb = 3.64 * 1000. # Blood specific heat [j/kg/k]
cair = 1.01 * 1000. # Air specific heat  [J./kg/K-]
emsk = 0.99 # Skin emissivity
emcl = 0.95 # Clothes emissivity
Lvap = 2.42 * 10. ** 6. # Latent heat of evaporation [J/Kg]
sigm = 5.67 * 10. ** (-8.) # Stefan-Boltzmann constant [W/(m2*K^(-4))]
eta = 0. # Body efficiency

# Initialisation of Temperature vector with respectively: Tcore, Tskin, Tcl
T = [38,40,40]
eps = 10**(-3) #Numerical precision

# Dichotomy browsning parameters
a = -40
b = 60
# Input data of the PET 
Ta=50 # Air temperature in oC
Tmrt=50 # Mean radiant temperature in oC
HR=50 # Air relative humidity %
v=0.1 # Wind velocity m/s
age = 35
sex = 1 # 1 for men and 2 for women
pos = 1
mbody = 75 #[Kg]
ht = 1.80 #[m]
p = 1013.25 #[hPa]
M_activity = 80 #[W] Activity
icl = 0.9# [clo] Clothing level

#initialisation pour le premier calcul
Tstable = resolution(Ta,Tmrt,HR,v,age,sex,ht,mbody,pos,M_activity,icl,T)[0]
print(PET(age, sex, ht, mbody, pos, M_activity, icl, Tstable, -30, 90, eps))