1 | c.......................................................................
   2 | c
   3 | c     calcul de la fonction "correction a la fonctionnelle poids"
   4 | c     pour tenir compte des homogeneites sur le trajet optique
   5 | c     modifiant le parametre de bande PHI (Curtis-Godson part II)
   6 | c     On calcul pdf(CG)/pdf(utilise)= fss(a,phicg)/fss(a/phi)
   7 | c     Il est deja tenu compte des modif de kgbar par kgskgbar(part I)
   8 | c
   9 | c     en entree : * a: k sur kgbar a l'emission (sans unite)
  10 | c                 * phi: parametre de forme a l'emission (")
  11 | c                 * phicg: idem mais de Curtis-Godson
  12 | c                 * tauba: la somme des kbar * l
  13 | c                 * kbarin: kbar du gaz dans maille d'emission ou sigma1
  14 | c                 * sigma1: longueur chemin dans maille emission
  15 | c
  16 | c                 *kbariin: kbar du gaz dans maille de reception
  17 | c                 *sigma_2: longueur chemin ds maille reception (inutile)
  18 | c                 *phi2: parametre de forme a la reception
  19 | c                 (les "deux" var ci-dessus servent pour CG en vol/vol)
  20 | c
  21 | c     en sortie : * la fonction correction pour un gaz
  22 | c                   a multiplier au poids (correction pdf * alpha)
  23 | c......................................................................
  24 | c     2000 juillet 9: Richard et Amaury
  25 | c     2000 octobre 12: " et " des modifs sur alpha
  26 | c.......................................................................
  27 |       double precision function cg(n1,n2,a,phi,
  28 |      & phicg,tauba,kbarin,sigma1,kbariin,sigma2,phi2)
  29 | !!     & phieq,ngazmx,nbandemx)
  30 | c.......................................................................
  31 | c      implicit double precision (a-h,o-z)
  32 |       implicit none
  33 |       include 'cecile.inc'  ! pout passer la varibale hetero
  34 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!      include 'propradia.inc'  !juste pour debug
  35 | c.......................................................................
  36 |       integer n1, n2 ,ngazmx,nbandemx
  37 |       double precision a, phi, phicg,kbarin,sigma1
  38 |       double precision kbariin, sigma2, phi2  !utilise pour vol/vol seul
  39 |       double precision toto
  40 |       double precision tauba
  41 |       double precision alpha, beta, df6f, dphicg6dl
  42 |       double precision phieqtronque  !sans le sigma_1
  43 |       double precision pdf1, pdf2
  44 | !      double precision phieq
  45 | !      dimension phieq(ngazmx, nbandemx)
  46 | c.......................................................................
  47 |       if (phi .eq. 0.d+0) then
  48 |          write (*,*) 'ERROR: division par phi nul (dans cg)'
  49 |          write (*,*) 'STOP'
  50 |          STOP
  51 | 
  52 |          elseif (a .eq. 0.d+0) then
  53 |          cg = 0.d+0
  54 |           !cg = 1.d+0
  55 |           !car la fonction qu'on inetgre a d'autre occasion
  56 |           !de s'annuler si elle doit l'etre
  57 |           !et ici ca fait une / 0 ensuite
  58 |          elseif (tauba .eq. 0.d+0) then
  59 |          cg = 1.d+0
  60 | 
  61 |          else
  62 |             !FLAG9
  63 |          toto = dexp(-(phicg-phi)*(a-1)**2/(2*a))
  64 |          cg = dsqrt(phicg/phi)*toto
  65 |       endif
  66 | !!      fss(phicg)/fss(phi)
  67 | ccccccc      call f_ss(phi,kbariin,1.d+0,a*kbariin,pdf1)
  68 | ccccccccc      call f_ss(phicg,kbariin,1.d+0,a*kbariin,pdf2)c
  69 | ccccccccc      cg = pdf2 /pdf1
  70 | cccccccccccc      if (pdf1.eq. 0.d+0) cg= 1.d+0
  71 |       !jusque la correction au pdf == hypothese de separabilite
  72 |       ! pour annuler cette correction on peut utiliser la ligne suivante
  73 |       !cg = 1.d+0
  74 |       alpha = 1.d+0 !pour l'instant
  75 | c..........................................................................
  76 |       if (hetero .eq. 2) then
  77 |          ! on y va pour curtis-godson
  78 |          ! alpha n'est pas egal a un 
  79 |          !rapport pdf * alpha
  80 |          !rapport pdf * (1-beta)
  81 |          !ATTENTION TROIS BRANCHES
  82 |          if (((n1.eq.0).or.(n1.eq.(n+1))).and.((n2.eq.0).or.
  83 |      &      (n2.eq.(n+1)))) then
  84 |             !la cas paroi/paroi
  85 |             !rien a faire cg*alpha avec alpha qui vaut un
  86 |             alpha = 1.d+0
  87 |             
  88 | 
  89 |             elseif ((n1.eq.0).or.(n1.eq.(n+1)).or.(n2.eq.0).or.
  90 |      &      (n2.eq.(n+1))) then
  91 |             !volume-paroi
  92 | c         if ((n1.eq.0).or.(n1.eq.(n+1)).or.(n2.eq.0).or.
  93 | c     &      (n2.eq.(n+1))) then
  94 |          ! pour (gaz paroi) ou (paroi paroi)
  95 |          !volume elementaire VERS paroi
  96 |          phieqtronque = phicg + (kbarin * sigma1 *(phicg - phi))/
  97 |      & (tauba-(kbarin * sigma1))
  98 | !!         if (n1.eq.2)  print *, 'phieqtronque, tauba, kbarin,sigma1'
  99 | !!         if (n1.eq.2) print *, phieqtronque, tauba, kbarin,sigma1
 100 |          !alpha, beta, df6f, dphicg6dl
 101 |         dphicg6dl = (kbarin *(tauba-(kbarin*sigma1))*(phi-phieqtronque))
 102 |      & /(tauba**2)
 103 |          !!Attention
 104 |          !!ceci pour maille a cote d'un paroi
 105 |          if (tauba .eq. (kbarin*sigma1)) dphicg6dl=0.d+0
 106 | !!       if (n1.eq.2)  print *, dphicg6dl
 107 |          df6f= 0.5d+0*((1/phicg)-((a-1)**2/a))
 108 | !!       if (n1.eq.2)  print *, df6f
 109 |          beta = df6f*(1/(kbarin*a))*dphicg6dl
 110 |          if (kbarin .eq. 0.d+0)  beta = 0.d+0
 111 |          !on n'emet pas de ray donc pas de correction
 112 |          alpha = (1 - beta)
 113 | !!!!!!!!!!!!!!!!!!!!         if (igaz.eq.3)  write(21,*) 'Moi  ', ibande,'      '
 114 | !!!!!!!!!!!!!!!!!!!     & ,  n1, n2, alpha !enlever proporadia
 115 | !         if ((n1.eq.0)) then !!!!!!!.and.(n2.eq.123)) then 
 116 |             !print *, 'phieqtronque,phicg,kbarin,sigma1,tauba'
 117 |             !print *, phieqtronque,phicg,kbarin,sigma1,tauba
 118 | !            print *, beta,df6f,dphicg6dl,kbarin,a,tauba
 119 | !            read *
 120 | !         endif
 121 | !!!         cg = cg * alpha
 122 | 
 123 |          !volume volume
 124 |          else
 125 |             alpha = 1.d+0
 126 |             !---MIEUX ci-dessous
 127 |             if (kbarin.eq.0.d+0) then 
 128 |                kbarin=1.d-20
 129 |             !!   phi= phieq(igaz,ibande)
 130 |             endif
 131 |             if (kbariin.eq.0.d+0) kbariin=1.d-20
 132 | !!!!            call cg_k_gg(a,phicg,tauba,phi,kbarin,phi2,kbariin,alpha)
 133 | !!            cg = cg * alpha
 134 |          endif
 135 | 
 136 | !!         alpha = 1.d+0 - ((phi-phicg)/(2*a*tauba))*
 137 | !!     & (1.d+0/phi - (a-1)**2.d+0/a)
 138 |          !print *, 'cg, alpha:  ',cg, alpha,'---',(phi-phicg),a,tauba,phi
 139 | !!         if (tauba .lt. 1.d-10) alpha = 1.d+0
 140 |          !print *, 'rrrrrrrrrrrrrr', cg, alpha
 141 |          !alpha = 1.d+0
 142 |          !print *, '>>>>>>>>>>>>>>>>>>CORRECTION:  ', cg, alpha, cg*alpha
 143 |          cg.f.html">cg = cg * alpha
 144 |          endif
 145 | 
 146 | c.......................................................................
 147 |       return
 148 |       end
 149 | c.......................................................................
 150 | 
 151 | 
 152 | 
 153 |