Field projections

A python routine to project two fields onto a regular 2D grid, according to a simplified adaptation of the algorithm of Aster's closest neighbours, also detailed on our page . The code is set in number of points around and can be executed in parallel, for example by using python multiprocessing.

##################################################################
# _______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
##################################################################
# ce sctript permet la projection de deux champs sur une grille  #
# régulière                                                      #
##################################################################
# Last modification :  14/08/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]com                                
##################################################################


import numpy as np
import math

def calc_parallel(argu_heure):

	print "quelle heure", argu_heure
	# fichier source
	dossier_src="./maillage/"
	nom_fichier=dossier_src+"coordonnees.txt"

	# ouverture
	coord=open(nom_fichier,"r")
	lignes=coord.readlines()
	coord.close()
	# on les met dans trois vecteurs uniques
	x=[]
	y=[]

	# comme les x,y sont en colonnes on les separe
	for line in lignes:
		columns = line.split(",") # separateur = virgule+espace(s)
		if not columns[0].startswith("#"): # verif inutile mais c'est cadeau
			x.append(float(columns[0]))
			y.append(float(columns[1]))

	# le type de flukx
	typeflux="diffus"
	dossier_src="./flux_"+typeflux+"/"
	dossier_dst="./flux_"+typeflux+"_projete/"
	suffixe="simulation_annuelle_"+typeflux+"_"
			
	# on definit les min / max des coordonnees
	xmin = min(x)
	xmax = max(x)
	ymin = min(y)
	ymax = max(y)
	#zmin = min(z)
	#zmax = max(z)

	#nb reel d'elements
	N_elts=len(x)
	#dx_reel_moyen=(xmax-xmin)/float(Nx_reel)
	
	# nb d'elements pour la grille reguliere
	Nx=180
	Ny=180

	# les pas reguliers d'espace (dx,dy)
	dx=(xmax-xmin)/float(Nx-1)
	dy=(ymax-ymin)/float(Ny-1)

	# print pour pouvoir en faire qque chose
	# print "Les coordonnees"
	# print "	xmin, xmax", xmin, xmax
	# print "	ymin, ymax", ymin, ymax
	# print "	dx,dy", dx, dy

	# les matrices contenant x,y reguliers (on peut faire sans en recodant)
	global grille_x
	global grille_y
	grille_x=np.ones((Nx,Ny))
	grille_y=np.ones((Nx,Ny))
	# une matrice pour verifier qu'on a bien rempli la grille
	global grille_verif
	grille_verif=np.zeros((Nx,Ny))

	# je mets ici l'exposant pour le lissage
	beta=0.333 # 2*0.75
	# le nb de points voisins pour voir
	nb_pts=4
	# une fc pour calculer la distance entre deux points 
	# (ca soulage les yeux dans la partie d'apres!)
	# il n'y a que 3 coordonnees car x et y ont le meme indice "i"
	def dist(i,x2,y2):
		dist_xy=math.sqrt(math.pow(x[i]-grille_x[x2][y2],2)+math.pow(y[i]-grille_y[x2][y2],2) )
		return dist_xy

	# creation grille reguliere
	for k in range(Nx):
		for l in range(Ny):
			grille_x[k][l] = k*dx+xmin
			grille_y[k][l] = l*dy+ymin
			
	for heure in range(argu_heure,argu_heure+1):
		# ponderations
		ponderation=np.zeros((Nx,Ny))
		somme_dist=np.zeros((Nx,Ny))
		champ=np.ones((Nx,Ny))*999.99
		
		print " le fichier", suffixe+str(heure)+".txt"
		
		nom_fichier_source=dossier_src+suffixe+str(heure)+".txt"
		nom_fichier_destination=dossier_dst+suffixe+str(heure)+".txt"
		fichier = open(nom_fichier_source, "r")
		z=map(float,fichier) #pour float > string
		fichier.close()
		# si z ne contient rien (pas de flux, genre la nuit)
		if z==[]:
			print "	(est vide)"
			fichier=open(nom_fichier_destination,"w")
			fichier.close()
		else:
			print " (n'est pas vide)"
			# balayer les points du nuage
			for i in range(N_elts):		
				# on trouve le point du maillage regulier le plus proche 
				# 	>> sur l'axe x
				if (x[i]- xmin) % dx > dx/2.:
					i_grille=1+int(math.floor( (x[i]-xmin)/dx ))
				else:
					i_grille=int(math.floor( (x[i]-xmin)/dx ))
				# 	>> sur l'axe y
				if (y[i]- ymin) % dy > dy/2.:
					j_grille=1+int(math.floor( (y[i]-ymin)/dy ))
				else:
					j_grille=int(math.floor( (y[i]-ymin)/dy ))

				# on stocke les valeurs de la grille qui ont ete remplies
				grille_verif[i_grille][j_grille]=1. # verif que touchee

				for m in range(-nb_pts,nb_pts+1):
			        for n in range(-nb_pts,nb_pts+1):
				        if i_grille+m>0 and j_grille+n>0 and j_grille+n<Ny and i_grille+m<Nx and grille_touched[i_grille+m][j_grille+n]==1:
        					distance=dist(i,i_grille+m,j_grille+n)
        					weight=math.exp(-math.pow(distance/dref,beta))
        					if beta==0:weight=1
        					ponderation[i_grille+m][j_grille+n]=ponderation[i_grille+m][j_grille+n]+weight*grille_touched[i_grille+m][j_grille+n]*z[i]
        					somme_dist[i_grille+m][j_grille+n]=  somme_dist[i_grille+m][j_grille+n]+weight*grille_touched[i_grille+m][j_grille+n]
        					grille_verif[i_grille+m][j_grille+n]=1

			print "Stockage dans la matrice"
			for i in range(Nx):
				for j in range(Ny):
					if grille_verif[i][j]==0: # si aucun point n'a ete touche
						zob=1
						#champ[i][j]=-5555.
					else:
						champ[i][j]= ponderation[i][j]/somme_dist[i][j]
			fichier=open(dossier_dst+suffixe+ str(heure)+".txt","w") 
			for i in range(Nx):
				for j in range(Ny):
					chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(champ[i][j])
					fichier.write(chaine+"\n") # retour ligne
				fichier.write("\n") # saut de ligne entre les x
			fichier.close()<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

A second, more advanced version consists of solving a minimization problem to project the field onto a grid (here regular) from the values of the point cloud. The code proposed below uses the library scipy.optimize.

# _______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
##################################################################
# ce sctript permet la projection de deux champs sur une grille  #
# régulière                                                      #
##################################################################
# Last modification :  16/01/2019                                #
##################################################################
# Copyright (C) 2018 Edouard WALTHER   / Antoine HUBERT          #  
# 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 time
import math
import csv
from scipy.optimize import minimize
debut=time.time()
alsace=True

algo="CG"
algo="SLSQP"
direction=1
algo="POWELL"
# nb de points en x et y
Nx=30
Ny=30
# la grille touchee (une matrice Nx*Ny remplie de 1 fait l'affaire si pas de "trous" dans votre maillage)
grille_touched=np.load("grille_touched_"+str(Nx)+ '.npy')
les_min_max=np.load("min_max"+ '.npy') #les min max des coordonnees
minxparoi1=les_min_max[0]
maxxparoi1=les_min_max[1]
minyparoi1=les_min_max[2]
maxyparoi1=les_min_max[3]

# quelques defs initiales pour les dossiers source et destination
if alsace==True:prefixe_dst="./vitesse_existant_hallHA"
else:prefixe_dst="./vitesse_existant_hallSM"

argu_heure=0
dossier_src="../../RESULTATS/hall_alsace/CSV_point/"
dossier_dst="./minimisation_matrice_beta_dref_flux_"+algo+"/"
#creation
if not os.path.isdir(dossier_dst):
	os.mkdir(dossier_dst)
#lecture des fichiers
if alsace==True:
	fichier_source=dossier_src+'flux_point.'+str(argu_heure)+".csv"

# def des coordonnees du champ : on prend un peu de marge
xmin = minxparoi1-4
xmax = maxxparoi1+4
ymin = minyparoi1-4
ymax = maxyparoi1+4
dx=(xmax-xmin)/float(Nx-1)
dy=(ymax-ymin)/float(Ny-1)
grille_x=np.ones((Nx,Ny))
grille_y=np.ones((Nx,Ny))

#une fc de calcul de la distance qui soulage les yeux
def dist(i,x2,y2):
	return math.sqrt(math.pow(x[i]-grille_x[x2][y2],2)+math.pow(y[i]-grille_y[x2][y2],2) )

# creation grille reguliere
for k in range(Nx):
	for l in range(Ny):
		grille_x[k][l] = k*dx+xmin
		grille_y[k][l] = l*dy+ymin
# lecture du fichiers
x,y,z = [],[],[]
with open(fichier_source,'r') as csvfile:
		reader = csv.reader(csvfile,delimiter=',',quoting = csv.QUOTE_NONNUMERIC)
		next(reader)
		for line in reader:
			x.append(line[1])
			y.append(line[2])
			z.append(line[0])

N_elts=len(x)
#-----------------------------------------------------
#
#	la fc de projection selon les plus proches voisins comme defini dans code Aster
#
def fc_projection_deg_1(vecteur_params):
	
	a=vecteur_params[0:Nx*Ny]
	b=vecteur_params[Nx*Ny:2*Nx*Ny]
	c=vecteur_params[2*Nx*Ny:3*Nx*Ny]
	beta=1.5 # exposant de "lissage" de la ponderation
	numerateur_dref=3.# le numerateur de dref=dx/numerateur_dref
	nb_pts=4 # cbien de points autour va-t-on chercher
	# ponderations
	dref=dx/numerateur_dref
	# balayer les points du nuage -----------------
	somme_moindre_ca=0.
	for i in range(N_elts):
	    # le carre le plus proche
		i_grille=int(math.floor((x[i]-xmin)/dx))
		j_grille=int(math.floor((y[i]-ymin)/dy))
		# on regarde les points autour
		for m in range(-nb_pts,nb_pts+1):
			for n in range(-nb_pts,nb_pts+1):
				if i_grille+m>0 and j_grille+n>0 and j_grille+n<Ny and i_grille+m<Nx and grille_touched[i_grille+m][j_grille+n]==1:
					distance=dist(i,i_grille+m,j_grille+n)
					weight=math.exp(-math.pow(distance/dref,beta))
					b_x = b[(i_grille+m)*Ny+j_grille+n]*grille_x[i_grille+m][j_grille+n]
					c_y = c[(i_grille+m)*Ny+j_grille+n]*grille_y[i_grille+m][j_grille+n]
					F= a[(i_grille+m)*Ny+j_grille+n] + b_x + c_y
					somme_moindre_ca=somme_moindre_ca+ weight*(F-z[i])**2
	return somme_moindre_ca

# les defs du 
a_tot=np.ones((3*Nx*Ny+2))
a_tot[-2]=1.5 #beta
a_tot[-1]=3. #numerateur_dref

# si on a une estimation precedente on peut la charger
#vec_precedent="./minimisation_matrice_beta_dref_flux_POWELL/vecteur_abc.npy"
#a_tot=np.load(vec_precedent)

# def des limites pour pour les coeffs a_k,b_k,c_k
#	et pour beta, num_dref
bnds=()
for k in range(3*Nx*Ny):
	bnds=bnds+((None,None),)
#limites pour beta 
bnds=bnds+((1,5),)
# ...et pour numerateur
bnds=bnds+((1,10),)

print("Debut minimisation pour l'algo " +algo + "...")

if algo=="POWELL":
	res=minimize(fc_projection_deg_1, a_tot, method='Powell', tol=1e-5, options={'disp':True, 'xtol':0.0001,'ftol':0.0001,'maxiter':10,'maxfev':30000})
if algo=="SLSQP":
	res=minimize(fc_projection_deg_1, a_tot, method='SLSQP', tol=1e-6, bounds=bnds, options={'disp':True,'ftol':0.00001,'maxiter':90})
if algo=="CG":
	res=minimize(fc_projection_deg_1, a_tot, method='CG', tol=1e-6, bounds=bnds, options={'disp':True, 'gtol':0.0001,'maxiter':10})
print(res.message)
print("Vecteur x", res.x)

#on sauvegarde le resultat
matrice=dossier_dst+"vecteur_abc.npy"
np.save(matrice, res.x)

#ensuite on ecrit sous le bon format pour gnuplot
vecteur_params=np.load(matrice)
a=vecteur_params[0:Nx*Ny]
b=vecteur_params[Nx*Ny:2*Nx*Ny]
c=vecteur_params[2*Nx*Ny:3*Nx*Ny]
print ("Stockage...")
fichier=open(dossier_dst+'flux_'+str(argu_heure)+".txt","w") 
for i in range(Nx):
	for j in range(Ny):
		# point touche ET dans la zone d'interet originelle
		if grille_touched[i][j]==0:
			chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(999.99)
		else:
			b_x = b[i*Ny+j]*grille_x[i][j]
			c_y = c[i*Ny+j]*grille_y[i][j]
			F= a[i*Ny+j] + b_x + c_y
			chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(F)
		fichier.write(str(chaine)+"\n") # retour ligne
	fichier.write("\n") # saut de ligne entre les x
fichier.close()
fin=time.time()

# calcul des temps passes
deltaT=round(fin-debut,1)
suffixe=" (s)"
if deltaT > 60:
	deltaT=round(deltaT/60,1)
	suffixe=" (min)"
if deltaT > 60:
	deltaT=round(deltaT/60,1)
	suffixe=" (h)"
print ("le calcul a duré "+str(deltaT)+suffixe +" pour l'algorithme " + algo)

Calculation of radiation view factors in urban areas

A Python tool developed during E. Medioni Master's internship (ENS Cachan) at AREP, based upon Radiance and gmsh.

These programs are used to determine the view factors of a STL geometry, obtained e.g. from Sketchup.

Download here (just add the gmsh executable in the "Resources" folder)

  • Mean Radiant Temperature

Calculation of the incident solar flux on a cylinder

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

S.E.T. calculation in Python

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

The corrected P.E.T.

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))