! ÉCOLE THÉMATIQUE INCENDIE – 30 MAI AU 04 JUIN 2015
! 
!   CALCUL DU RAYONNEMENT ÉMIS PAR UNE FLAMME CYLINDRIQUE VERTICALE ET REÇU PAR LE SOL
!   PAR LA MÉTHODE DE MONTE CARLO
! CONTACTS:     MATTHIEU.DEGENNARO@NOVELTIS.FR
!               MAXIME.MENSE@GMAIL.COM
!               BERNARD.PORTERIE@UNIV-AMU.FR
! AIX-MARSEILLE UNIVERSITÉ, CNRS, IUSTI UMR 7343, 13453 MARSEILLE, FRANCE
!
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,K,L,N,M,II,JJ,IMAX,JMAX,NQUANTA
PARAMETER (IMAX=101,JMAX=101)	! NOMBRE DE VOLUMES DE CONTRÔLE DU DOMAINE AU SOL
PARAMETER (M=41,L=41)	    ! MAILLAGE FLAMME
REAL*8 LX,LY,HF,DYP,DZF,YP,ZP,TETA,FI,RAN,PI,RC,SIGMA,  &
       Q0,FLUX(IMAX,JMAX),XP,DXP,R,ZARR,YARR,XARR,	&
	   RF,DS,XF(M,L),YF(M,L),ZF(M,L)
REAL*8 UX,UY,UZ,VX,VY,VZ,NORMX,NORMY,NORMZ,XC,YC,ZC,	&
	   NORME1,NORME2,NORME,X1,Y1,Z1,X3,Y3,Z3,DGAMA,GAMA,PE
!===========================================================================
HF=0.5D0 ! HAUTEUR DE FLAMME
RF=0.1D0 ! RAYON DE LA FLAMME
LX=2.D0 ! LONGUEUR DU DOMAINE AU SOL
LY=2.D0 ! LARGEUR DU DOMAINE AU SOL
PE=100.D0 ! POUVOIR ÉMISSIF (FACTEUR DE FORME: METTRE PE=1.d0)
NQUANTA=1.D8 ! NOMBRE DE QUANTA PAR MÈTRE CARRÉ DE FLAMME
Q0=PE/NQUANTA !PUISSANCE TRANSPORTÉE PAR CHAQUE QUANTUM
R=100.D0 ! RAYON DE LA SPHÈRE ENVELOPPE (ARBITRAIRE, R>>LX)
!===========================================================================
SIGMA=5.67D-8
PI=DACOS(-1.D0)
FLUX=0.D0
!===========================================================================
! INITIALISATION DU GÉNÉRATEUR DE NOMBRES ALÉATOIRES
!===========================================================================
CALL RANDOM_SEED
!===========================================================================
! MAILLAGE DU DOMAINE AU SOL 
!===========================================================================
DXP=LX/IMAX
DYP=LY/JMAX
!===========================================================================
! MAILLAGE DE LA FLAMME 
!===========================================================================
DGAMA=2.D0*PI/(M-1)
DZF=HF/(L-1)
DO I=1,M
GAMA=(I-1)*DGAMA
RC=RF
DO J=1,L
XF(I,J)=RC*DCOS(GAMA)
YF(I,J)=RC*DSIN(GAMA)
ZF(I,J)=(J-1)*DZF
END DO
END DO
!---ECRITURE DE LA FLAMME DANS FICHIER
OPEN (UNIT=18,FILE='FLAMME.DAT')
WRITE(18,*)'TITLE="3D"'
WRITE(18,*)'VARIABLES="X","Y","Z"'
WRITE(18,*)'ZONE I=',L,'J=',M
WRITE(18,2003)((XF(I,J),YF(I,J),ZF(I,J),J=1,L),I=1,M)
!===========================================================================
! CALCUL DE LA SURFACE ÉLÉMENTAIRE ET DU NOMBRE DE QUANTA ÉMIS PAR CETTE SURFACE
!===========================================================================
DS=RF*DGAMA*DZF
N=INT(NQUANTA*DS)
!===========================================================================
! BOUCLE SUR LES ÉLÉMENTS D’ÉMISSION
!===========================================================================
DO I=1,M-1
PRINT *,'I=',I
DO J=1,L-1
        ! CALCUL DES COORDONNÉES DE LA NORMALE (NORMX, NORMY, NORMZ) ET DU 
        ! BARYCENTRE (XC, YC, ZC) DE CHAQUE ÉLÉMENT DE SURFACE DE FLAMME
        UX=XF(I+1,J)-XF(I,J)
		UY=YF(I+1,J)-YF(I,J)
		UZ=ZF(I+1,J)-ZF(I,J)
		NORME1=DSQRT(UX*UX+UY*UY+UZ*UZ)
		UX=UX/NORME1
		UY=UY/NORME1
		UZ=UZ/NORME1
		VX=XF(I,J+1)-XF(I,J)
		VY=YF(I,J+1)-YF(I,J)
		VZ=ZF(I,J+1)-ZF(I,J)
		NORME2=DSQRT(VX*VX+VY*VY+VZ*VZ)
		VX=VX/NORME2
		VY=VY/NORME2
		VZ=VZ/NORME2
		NORMX=UY*VZ-UZ*VY
		NORMY=UZ*VX-VZ*UX
		NORMZ=UX*VY-VX*UY
		NORME=DSQRT(NORMX**2+NORMY**2+NORMZ**2)
		NORMX=NORMX/NORME
		NORMY=NORMY/NORME
		NORMZ=NORMZ/NORME
            ! CALCUL DES COORDONNEES DU BARYCENTRE DE DS, POINT D'ÉMISSION DES QUANTA
		    XC=(XF(I+1,J+1)+XF(I,J+1)+XF(I,J)+XF(I+1,J))/4.D0
		    YC=(YF(I+1,J+1)+YF(I,J+1)+YF(I,J)+YF(I+1,J))/4.D0
		    ZC=(ZF(I+1,J+1)+ZF(I,J+1)+ZF(I,J)+ZF(I+1,J))/4.D0
                    !BOUCLE SUR LES QUANTA
		            PHOTON : DO K=1,N
            	        !CALCUL DES ANGLES D’ÉMISSION : THETA ET PHI
			        CALL RANDOM_NUMBER(RAN)
                    TETA=DASIN(DSQRT(RAN))
			        CALL RANDOM_NUMBER(RAN)
                    FI=2*PI*RAN
			        !CALCUL DES COORDONNÉES DU POINT D’INTERSECTION 
                    !DE LA DIRECTION D’ÉMISSION AVEC UNE SPHÈRE DE RAYON R
                    X3=R*DSIN(TETA)*DCOS(FI)
			        Y3=R*DSIN(TETA)*DSIN(FI)
			        Z3=R*DCOS(TETA)
			        !CHANGEMENT DE REPÈRE
                    X1=X3*UX+Y3*VX+Z3*NORMX
			        Y1=X3*UY+Y3*VY+Z3*NORMY
			        Z1=X3*UZ+Y3*VZ+Z3*NORMZ
                    !CALCUL DES COORDONNÉES DU POINT D'ARRIVÉE DANS LE REPÈRE FIXE
			        XARR=XC+X1
			        YARR=YC+Y1
			        ZARR=ZC+Z1
			        !	SI LE QUANTUM TOUCHE LE SOL (ZARR <= 0), ON CALCULE LES COORDONNÉES (XP,YP)
                    !DU POINT D’IMPACT, ET ON RAJOUTE Q0 AU PAVÉ DU SOL TOUCHÉ PAR LE QUANTUM
			        IF(ZARR.LE.0.D0)THEN
                        ZP=0.D0
				        XP=((ZP-ZC)/(ZARR-ZP)*XARR+XC)/(1+(ZP-ZC)/(ZARR-ZP))
				        YP=((ZP-ZC)/(ZARR-ZP)*YARR+YC)/(1+(ZP-ZC)/(ZARR-ZP))
				        IF(XP.GE.-LX/2.D0.AND.XP.LE.LX/2.D0.AND.YP.GE.-LY/2.D0.AND.YP.LE.LY/2.D0) THEN
				        II=INT((XP+LX/2.D0)/DXP)+1
				        JJ=INT((YP+LX/2.D0)/DYP)+1
				        FLUX(II,JJ)=FLUX(II,JJ)+Q0
				        ENDIF 
			        ENDIF
		            ENDDO PHOTON
ENDDO
ENDDO
!===========================================================================
! ENREGISTREMENT SUR FICHIER
!===========================================================================
OPEN (UNIT=17,FILE='FLUX_SOL.DAT')
WRITE(17,*)'TITLE="FLUX RAY AU SOL"'
WRITE(17,*)'VARIABLES="X","Y","FLUX"'
WRITE(17,*)'ZONE I=',IMAX,'J=',JMAX
WRITE(17,2003) ((-LX/2.D0+(2*I-1)*DXP/2.D0,-LY/2.D0+(2*J-1)*DYP/2.D0,FLUX(I,J)/(DXP*DYP),I=1,IMAX),J=1,JMAX)
!===========================================================================
2002 FORMAT(2(E11.4,1X))
2003 FORMAT(3(E11.4,1X))
2004 FORMAT(4(E11.4,1X))
!===========================================================================
STOP
END