! ÉCOLE THÉMATIQUE INCENDIE – 30 MAI AU 04 JUIN 2015
! 
!   CALCUL DU RAYONNEMENT ÉMIS PAR UN PANNEAU VERTICAL 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, PARAMETER :: NQUANTA=1.D6 ! NOMBRE DE QUANTA/M² DE PANNEAU
INTEGER, PARAMETER :: IMAX=100,JMAX=100	! NOMBRE DE VOLUMES DE CONTRÔLE DU DOMAINE AU SOL
INTEGER, PARAMETER :: JMAXP=40,KMAXP=40	! NOMBRE DE VOLUMES DE CONTRÔLE DU PANNEAU
INTEGER I,J,K,N,II,JJ,KK,NDSP
DOUBLE PRECISION LARGEUR_DOMAINE,LONGUEUR_DOMAINE,LARGEUR_PANNEAU,HAUTEUR_PANNEAU
DOUBLE PRECISION DX,DY,DYP,DZP,TETA,FI,RAN,PI,Q0,P(IMAX,JMAX),X(IMAX),Y(JMAX),YP(JMAXP),ZP(KMAXP)
DOUBLE PRECISION XC,YC,ZC,XARR,YARR,PE,R,DS,DSP
DOUBLE PRECISION EW_SOL
!===========================================================================
LONGUEUR_DOMAINE=20.D0  ! LONGUEUR DU DOMAINE AU SOL
LARGEUR_DOMAINE=20.D0   ! LARGEUR DU DOMAINE AU SOL
LARGEUR_PANNEAU=10.D0   ! LONGUEUR DU PANNEAU
HAUTEUR_PANNEAU=10.D0   ! HAUTEUR DU PANNEAU
EW_SOL=1.D0             ! ÉMISSIVITÉ DU SOL
PE=100.D0               ! POUVOIR ÉMISSIF EN KW/M²
Q0=PE/NQUANTA           ! PUISSANCE TRANSPORTÉE PAR UN QUANTUM
!===========================================================================
PI=DACOS(-1.D0)
!===========================================================================
CALL RANDOM_SEED
P=0.D0
!===================================================================
! MAILLAGES DU DOMAINE AU SOL ET DU PANNEAU RADIANT
!===================================================================
DX=LONGUEUR_DOMAINE/IMAX
DY=LARGEUR_DOMAINE/JMAX
DYP=LARGEUR_PANNEAU/JMAXP
DZP=HAUTEUR_PANNEAU/KMAXP
!--MAILLAGE DU DOMAINE AU SOL
DO II=1,IMAX    
	X(II)=(II-0.5D0)*DX
    !PRINT *,'I,X=',II,X(II)
ENDDO
DO JJ=1,JMAX
	Y(JJ)=-LARGEUR_DOMAINE/2.D0+(JJ-0.5D0)*DY
    !PRINT *,'JJ,Y=',JJ,Y(JJ)
ENDDO
!
DO JJ=1,JMAXP
	YP(JJ)=-LARGEUR_PANNEAU/2.D0+(JJ-0.5D0)*DYP
    !PRINT *,'JJ,YP=',JJ,YP(JJ)
ENDDO
DO KK=1,KMAXP
	ZP(KK)=(KK-0.5D0)*DZP
    !PRINT *,'KK,ZP=',KK,ZP(KK)
ENDDO
!===================================================================
DS=DX*DY    ! SURFACE ELEMENTAIRE DOMAINE SOL
DSP=DYP*DZP ! SURFACE ELEMENTAIRE PANNEAU
!===================================================================
! CALCUL DU NOMBRE DE QUANTA ÉMIS PAR ÉLÉMENT DE SURFACE DU PANNEAU 
!===================================================================
NDSP=INT(NQUANTA*DSP)
!===================================================================
! BOUCLE SUR LES ÉLÉMENTS DE SURFACE D’ÉMISSION DU PANNEAU
!===================================================================
DO J=1,JMAXP
PRINT *,'J=',J
DO K=1,KMAXP
! COORDONNÉES DU POINT D'ÉMISSION C  
XC=0.D0     
YC=YP(J)
ZC=ZP(K)
            ! BOUCLE SUR LES QUANTA ÉMIS PAR ÉLÉMENT DE SURFACE DU PANNEAU
            PHOTON : DO N=1,NDSP
            ! ANGLES D'ÉMISSION
            CALL RANDOM_NUMBER(RAN)
            TETA=DASIN(DSQRT(RAN))
            CALL RANDOM_NUMBER(RAN)
            FI=2*PI*RAN
            ! DISTANCE ENTRE LE POINT D'ÉMISSION ET LE SOL
            R=-ZC/(DSIN(TETA)*DSIN(FI))
            ! ON ÉLIMINE LA SOLUTION R<0
            IF(R.LT.0.D0) GOTO 120
            ! COORDONNÉES DU POINT D'IMPACT AU SOL
            XARR=XC+R*DCOS(TETA)
            YARR=YC+R*DSIN(TETA)*DCOS(FI)
            ! ON TESTE SI LE POINT D'IMPACT EST DANS LE DOMAINE AU SOL
            IF(XARR.GE.0.D0.AND.XARR.LE.LONGUEUR_DOMAINE.AND.   &
                    YARR.GE.-LARGEUR_DOMAINE/2.D0.AND.YARR.LE.LARGEUR_DOMAINE/2.D0) THEN
            ! RÉFLEXION OU ABSORPTION AU SOL ? 
            CALL RANDOM_NUMBER(RAN)
            IF(RAN.GE.EW_SOL) GOTO 120 
            ! ON RAJOUTE Q0 EN CAS D'ABSORPTION DU QUANTUM PAR LE PAVÉ II,JJ DU DOMAINE AU SOL
            II=INT(XARR/DX)+1
            JJ=INT((YARR+LARGEUR_DOMAINE/2.D0)/DY)+1
            P(II,JJ)=P(II,JJ)+Q0
            ENDIF 
            120             CONTINUE
            ENDDO PHOTON
ENDDO
ENDDO
!===================================================================
! ENREGISTREMENT SUR FICHIERS DU FLUX REÇU AU SOL ET LE LONG DE L'AXE X
!===================================================================
OPEN (UNIT=17,FILE='FLUX.DAT')
WRITE(17,*)'TITLE="FLUX AU SOL"'
WRITE(17,*)'VARIABLES="X","Y","FLUX(KW/M2)"'
WRITE(17,*)'ZONE I=',IMAX,'J=',JMAX
WRITE(17,2003) ((X(I),Y(J),P(I,J)/DS,I=1,IMAX),J=1,JMAX)
!
OPEN (UNIT=21,FILE='FLUX_VS_X.DAT')
WRITE(21,*)'TITLE="FLUX SELON X"'
WRITE(21,*)'VARIABLES="X","FLUX(KW/M2)"'
IF(MOD(JMAX,2).EQ.0) THEN
WRITE(21,2002) (X(I),(P(I,JMAX/2)+P(I,JMAX/2+1))/2.D0/DS,I=1,IMAX)
ELSE
WRITE(21,2002) (X(I),P(I,(JMAX+1)/2)/DS,I=1,IMAX)
ENDIF
!===================================================================
2002 FORMAT(2(E11.4,1X))
2003 FORMAT(3(E11.4,1X))
2004 FORMAT(4(E11.4,1X))
STOP
END