Sensibilités (programmes 8 et 8-abs)

Programme 8 Reprendre le programme 3 en évaluant également la sensibilité de $\left \langle L^2 \right \rangle$ à $k_d=\frac{1}{\lambda_d}$, ainsi que la sensibilité de $\left \langle 1 - exp\left(-\frac{L}{\lambda_a}\right) \right \rangle$ à $k_a=\frac{1}{\lambda_a}$.

Exemple de solution pour le carré des longueurs (programmé en fortran) :

c---------------------------------------------------------------------
c Calcul de la sensibilité à k de la moyenne des carrés de longueurs
c---------------------------------------------------------------------
      PROGRAM prog4
      IMPLICIT NONE
      INTEGER iseed
      COMMON/comsee/iseed
      DOUBLE PRECISION pi
      PARAMETER(pi=3.141592654)
c---------------------------------------------------------------------
      INTEGER itir,ntir
      DOUBLE PRECISION w,som_w,som_w2,moy_w,moy_w2,sig_w,sig_moy_w
      DOUBLE PRECISION cote,lambda
      DOUBLE PRECISION x,y,theta,L
      DOUBLE PRECISION r,xold,yold,Lold,d,dsecx,dsecy
      DOUBLE PRECISION k,k0,dk
      DOUBLE PRECISION som_lambda_moins_d
      INTEGER ik,nk
      DOUBLE PRECISION sensib_k_moy_L2,sig_sensib_k_moy_L2
c---------------------------------------------------------------------
      OPEN(10,FILE='sorties.out')
      PRINT *,'germe du generateur aleatoire : '
      READ *,iseed
c---------------------------------------------------------------------

      cote=0.3
      nk=5
      k0=0.
      dk=1./cote
      ntir=100000

      DO ik=1,nk
         k=k0+ik*dk
         lambda=1./k

         som_w=0.
         som_w2=0.
         DO itir=1,ntir

            L=0.
            som_lambda_moins_d=0.

            x=0.
            CALL rand_uniforme(r)
            y=r*cote
            CALL rand_uniforme(r)
            theta=asin(2.*r-1.)

 1          CONTINUE
            
            xold=x
            yold=y
            Lold=L
         
            CALL rand_uniforme(r)
            d=-lambda*log(1-r)
            x=x+d*cos(theta)
            y=y+d*sin(theta)
            L=L+d
            som_lambda_moins_d=som_lambda_moins_d+(lambda-d)

            IF ((x.GE.0).AND.(x.LE.cote)
     &           .AND.(y.GE.0).AND.(y.LE.cote)) THEN
               CALL rand_uniforme(r)
               theta=r*2*pi
               GOTO 1
            ELSE
               IF (cos(theta).GE.0.) THEN
                  dsecx=(cote-xold)/cos(theta)
               ELSE
                  dsecx=(0.-xold)/cos(theta)
               ENDIF
               IF (sin(theta).GE.0.) THEN
                  dsecy=(cote-yold)/sin(theta)
               ELSE
                  dsecy=(0.-yold)/sin(theta)
               ENDIF
               IF (dsecx.LE.dsecy) THEN
                  d=dsecx
               ELSE
                  d=dsecy
               ENDIF
               x=xold+d*cos(theta)
               y=yold+d*sin(theta)
               L=Lold+d
            ENDIF

            w=som_lambda_moins_d*L**2
            
            som_w=som_w+w
            som_w2=som_w2+w**2
         ENDDO
         moy_w=som_w/ntir
         moy_w2=som_w2/ntir
         sig_w=sqrt(moy_w2-moy_w**2)
         sig_moy_w=sig_w/sqrt(dble(ntir))

         sensib_k_moy_L2=moy_w
         sig_sensib_k_moy_L2=sig_moy_w
         WRITE(10,*) k,sensib_k_moy_L2,sig_sensib_k_moy_L2

      ENDDO

c---------------------------------------------------------------------
      CLOSE(10)
      STOP
      END
c---------------------------------------------------------------------

Exemple de solution pour l'absorption (programmé en fortran) :

c---------------------------------------------------------------------
c Calcul de la sensibilité au coefficient d'absorption (ka) de la
c fraction absorbée en fonction de l'inverse du libre parcours moyen
c de diffusion (k)
c---------------------------------------------------------------------
      PROGRAM prog3
      IMPLICIT NONE
      INTEGER iseed
      COMMON/comsee/iseed
      DOUBLE PRECISION pi
      PARAMETER(pi=3.141592654)
c---------------------------------------------------------------------
      INTEGER itir,ntir
      DOUBLE PRECISION w,som_w,som_w2,moy_w,moy_w2,sig_w,sig_moy_w
      DOUBLE PRECISION cote,lambda
      DOUBLE PRECISION x,y,theta,L
      DOUBLE PRECISION r,xold,yold,Lold,d,dsecx,dsecy
      DOUBLE PRECISION ka,k,k0,dk
      INTEGER ik,nk
      DOUBLE PRECISION sensib_ka_frac_abs,sig_sensib_ka_frac_abs
c---------------------------------------------------------------------
      OPEN(10,FILE='sorties.out')
      PRINT *,'germe du generateur aleatoire : '
      READ *,iseed
c---------------------------------------------------------------------

      cote=0.3
      ka=0.01/cote
      nk=5
      k0=0.
      dk=1./cote
      ntir=10000

      DO ik=1,nk
         k=k0+ik*dk
         lambda=1./k

         som_w=0.
         som_w2=0.
         DO itir=1,ntir

            L=0.

            x=0.
            CALL rand_uniforme(r)
            y=r*cote
            CALL rand_uniforme(r)
            theta=asin(2.*r-1.)

 1          CONTINUE
            
            xold=x
            yold=y
            Lold=L
         
            CALL rand_uniforme(r)
            d=-lambda*log(1-r)
            x=x+d*cos(theta)
            y=y+d*sin(theta)
            L=L+d

            IF ((x.GE.0).AND.(x.LE.cote)
     &           .AND.(y.GE.0).AND.(y.LE.cote)) THEN
               CALL rand_uniforme(r)
               theta=r*2*pi
               GOTO 1
            ELSE
               IF (cos(theta).GE.0.) THEN
                  dsecx=(cote-xold)/cos(theta)
               ELSE
                  dsecx=(0.-xold)/cos(theta)
               ENDIF
               IF (sin(theta).GE.0.) THEN
                  dsecy=(cote-yold)/sin(theta)
               ELSE
                  dsecy=(0.-yold)/sin(theta)
               ENDIF
               IF (dsecx.LE.dsecy) THEN
                  d=dsecx
               ELSE
                  d=dsecy
               ENDIF
               x=xold+d*cos(theta)
               y=yold+d*sin(theta)
               L=Lold+d
            ENDIF

            w=L*exp(-ka*L)
            
            som_w=som_w+w
            som_w2=som_w2+w**2
         ENDDO
         moy_w=som_w/ntir
         moy_w2=som_w2/ntir
         sig_w=sqrt(moy_w2-moy_w**2)
         sig_moy_w=sig_w/sqrt(dble(ntir))

         sensib_ka_frac_abs=moy_w
         sig_sensib_ka_frac_abs=sig_moy_w
         WRITE(10,*) k,sensib_ka_frac_abs,sig_sensib_ka_frac_abs

      ENDDO

c---------------------------------------------------------------------
      CLOSE(10)
      STOP
      END
c---------------------------------------------------------------------

Exemple d'extension en géométrie tridimensionnelle (sous EDStar) : en cours.