1 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   2 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   3 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   4 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   5 | c
   6 | c.......................................................................
   7 | c
   8 | c  inverse gaussian cumulative
   9 | c
  10 | c  in :  * $\phi$                ---> phgicd
  11 | c        * $\mu$                 ---> umgicd
  12 | c        * $x$                   ---> xgicd
  13 | c
  14 | c  out : * $\int_0^x \Im (t ; \mu , \mu \phi ) dt$
  15 | c                                ---> cdgicd
  16 | c
  17 | c.......................................................................
  18 |       subroutine gicd(phgicd,umgicd,xgicd,cdgicd)
  19 | c.......................................................................
  20 |       implicit double precision (a-h,o-z)
  21 | c.......................................................................
  22 |       !print *, 'GICD: ', phgicd,umgicd, xgicd
  23 |       gicd1=xgicd/umgicd
  24 |       gicd2=dsqrt(phgicd/gicd1)
  25 |       gicd3=-gicd2*(1.d0-gicd1)
  26 |       gicd4=-gicd2*(1.d0+gicd1)
  27 |       call stcd(gicd3,gicd5)
  28 |       call stcd(gicd4,gicd6)
  29 |       cdgicd=gicd5+dexp(2.d0*phgicd)*gicd6
  30 | 
  31 |       if (cdgicd .ge. 1.d+0) then
  32 |           cdgicd = 1.d+0 - 1.d-8
  33 |       endif
  34 |       !print *, 'toto:= ', gicd6, gicd5, phgicd, dexp(2.d0*phgicd)
  35 |       !print *, 'fim:= ', gicd5+dexp(2.d0*phgicd)*gicd6
  36 | c....tjours pour soigner toto
  37 |       if (phgicd .ge. 400) then
  38 |       !if (phgicd .ge. 600) then
  39 |          !phi_star est tres grand la proba de k = un dirac
  40 |          !la cumulative vaut alors
  41 |          if (xgicd .lt. umgicd) then
  42 |             !zero en dessous du dirac
  43 |             cdgicd = 0.d+0
  44 |             else
  45 |             !un au dessus du dirac
  46 |             cdgicd = 1.d+0
  47 |          endif
  48 |       endif
  49 |       !print *, 'GICD: ', cdgicd, dexp(2.d0*phgicd), phgicd
  50 | c.......................................................................
  51 |       return
  52 |       end