C[name] th36
C[parameter_names]  phi, n, nw, rhow, rho1, rho2, sig1, sig2, Rmain, zz, Rshift, bckgr
C[parameter_init_values] 0.01, 3, 0.1, -0.56e10, 2e10, 3e10, 50.0, 20.0, 800.0, 0.2, 10.0, 0.01
C[parameter_init_range] [0.0..1.0], [1..100], [0.0..1.0], [-1e10..10e10], [-1e10..10e10], [-1e10..10e10], [0.0..10000.0], [0.0..10000], [0.0..100000.0], [0.0..10000], [0.0..10000.0], [0.0..1.0] 
C[parameter_units] [1], [1], [1], [1/cm^2], [1/cm^2], [1/cm^2], [A], [A], [A], [?1], [A], [1/cm] 
C[parameter_info]  [Konzentration Polymer], [Anzahl der Schichten], [width of n-distribution], [SLD Wasser], [SLD der zentralen Schicht (1)], [SLD der Dekoration (2)], [Dicke der zentralen Schicht (1)], [Dicke der Dekoration (2)], [Main Radius], [Z - Schultz distribution], [Maximum possible shift], [background]
![info]  "Multilamellar vesicle model",, "Model author: Henrich Frielinghaus",, "https://journals.aps.org/pre/abstract/10.1103/PhysRevE.76.051603",, "Fortran code: Henrich Frielinghaus",, "qtikws adaptation: Vitaliy Pipich"

      function th36(x,xphi,xn,xnw,xrhow,xrho1,xrho2,xsig1,xsig2,
     *              xRmain,xzz,xRshift,xbckgr)
!
      implicit none
C
      double precision th36
      double precision x
C
      double precision Pi
      parameter (Pi = 3.14159265358979323846264338328d0)
C
      double precision xphi        ! Konzentration Polymer
      double precision xn          ! Anzahl der Schichten
      double precision xnw         ! width of n-distribution
      double precision xrhoW       ! SLD   Wasser
      double precision xrho1,xrho2 ! SLD   der zentralen Schicht (1) und Dekoration (2)
      double precision xsig1,xsig2 ! Dicke der zentralen Schicht (1) und Dekoration (2)
      double precision xRmain      ! Main Radius
      double precision xzz         ! Z - Schultz distribution
      double precision xRshift     ! Maximum possible shift
      double precision xbckgr      ! background
c
      double precision q           ! Q-Vektor
      double precision phi         ! Konzentration Polymer
      double precision n           ! Anzahl der Schichten
      double precision nw          ! width of n-distribution
      double precision rhoW        ! SLD   Wasser
      double precision rho1,rho2   ! SLD   der zentralen Schicht (1) und Dekoration (2)
      double precision sig1,sig2   ! Dicke der zentralen Schicht (1) und Dekoration (2)
      double precision Rmain       ! Main Radius
      double precision zz          ! Z - Schultz distribution
      double precision Rshift      ! Maximum possible shift
      double precision bckgr       ! background
C
      double precision drho1,drho2 ! Kontrast der zentralen Schicht (1) und Dekoration (2) zum L"osungsmittel
C
      double precision bmin,bmax   ! Integrationsgrenzen
      double precision eps         ! Fehler
      integer          maxit       ! Maximale Iterationen
c     double precision erracc      ! Tats"achlicher Fehler
      double precision aint        ! Integral
      integer          i
      double precision R0          ! actual value
C
      common /VesPar4/q,phi,n,nw,drho1,drho2,sig1,sig2,Rmain,zz,Rshift,
     x                bckgr
C
      double precision inVesicle4a
      external inVesicle4a
C
C     initialisation of names and number of parameters
C
        q     = x
        phi   = xphi
        n     = xn
        nw    = xnw
        rhow  = xrhow
        rho1  = xrho1
        rho2  = xrho2
        sig1  = xsig1
        sig2  = xsig2
        Rmain = xRmain
        zz    = xzz
        Rshift= xRshift
        bckgr = xbckgr
C
        n = max(1d0,n)
C
        drho1  = rho1 - rhoW
        drho2  = rho2 - rhoW
C
        if (zz.le.-.99d0) then

        th36 = real(inVesicle4a(Rmain))

        else

        if (abs(Rshift).lt.0.5d0*abs(Rmain)) Rshift = 0.5d0*abs(Rmain)   !! practical limit !?...

        aint  = 0d0

        bmin  = Rmain / sqrt(zz+1d0) / 3d0
        bmax  = Rmain / sqrt(zz+1d0)
        maxit = 100
        eps   = (bmax-bmin)/maxit

        do i=1,maxit
         R0   = (i-0.5d0)*eps + bmin
         aint = aint + inVesicle4a(R0)*eps
        enddo

        bmin  = Rmain * sqrt(zz+1d0)
        bmax  = Rmain * sqrt(zz+1d0) * 3d0
        maxit = 100
        eps   = (bmax-bmin)/maxit

        do i=1,maxit
         R0   = (i-0.5d0)*eps + bmin
         aint = aint + inVesicle4a(R0)*eps
        enddo

        bmin  = Rmain / sqrt(zz+1d0)
        bmax  = Rmain * sqrt(zz+1d0)
        maxit = 1000
        eps   = (bmax-bmin)/maxit

        do i=1,maxit
         R0   = (i-0.5d0)*eps + bmin
         aint = aint + inVesicle4a(R0)*eps
        enddo

        th36 = real(aint)

        endif
C
      return
      end
C
C-------------------------------------------------------------
C
      double precision function inVesicle4a(R0)

      implicit none

      double precision R0         ! Radius einer Schicht (verteilt)

      double precision q           ! Q-Vektor
      double precision phi         ! Konzentration Polymer
      double precision n           ! Anzahl der Schichten
      double precision nw          ! width of n-distribution
      double precision drho1,drho2 ! Kontrast der zentralen Schicht (1) und Dekoration (2) zum L"osungsmittel
      double precision sig1,sig2   ! Dicke    der zentralen Schicht (1) und Dekoration (2)
      double precision Rmain       ! Main Radius
      double precision zz          ! Z - Schultz distribution
      double precision Rshift      ! Maximum possible shift
      double precision bckgr       ! background

      integer          n1,n2,n3    ! n-variables
      double precision wt,sw       ! weights

      double precision inVesicle4
      external inVesicle4

      common /VesPar4/q,phi,n,nw,drho1,drho2,sig1,sig2,Rmain,zz,Rshift,
     x                bckgr

C--------------------------------------------

      if (nw.lt.0.5d0) then

        n1 = int(n)
        inVesicle4a = (n1+1d0-n)*inVesicle4(R0,dble(n1))
     x              + (n-n1)    *inVesicle4(R0,dble(n1)+1d0)

      else

        n1 = max(int(n-2d0*nw),1)
        n2 =     int(n+2d0*nw+1d0)

        inVesicle4a = 0d0
        sw = 0d0

        do n3=n1,n2
          wt = exp(-0.5d0*(n3-n)*(n3-n)/(nw*nw))
          inVesicle4a = inVesicle4a + wt*inVesicle4(R0,dble(n3))
          sw = sw + wt
        enddo

        inVesicle4a = inVesicle4a/sw

      endif

      end


C
C-------------------------------------------------------------
C
      double precision function inVesicle4(R0,n)

      implicit none

      double precision R0         ! Radius einer Schicht (verteilt)
      double precision n          ! number of layers

      double precision Pi
      parameter (Pi = 3.14159265358979323846264338328d0)

      double precision eps
      parameter (eps = 1d-6)

      double precision Vorf        ! Vorfaktor
      double precision Fv          ! Formfaktor / sqrt(2Pi)
      double precision Sv          ! Strukturfaktor (des geschichteten Vesikels)
      double precision f_r         ! Schultz distribution

      double precision q           ! Q-Vektor
      double precision phi         ! Konzentration Polymer
      double precision no          ! Anzahl der Schichten
      double precision nw          ! width of n-distribution
      double precision drho1,drho2 ! Kontrast der zentralen Schicht (1) und Dekoration (2) zum L"osungsmittel
      double precision sig1,sig2   ! Dicke    der zentralen Schicht (1) und Dekoration (2)
      double precision Rmain       ! Main Radius
      double precision zz          ! Z - Schultz distribution
      double precision Rshift      ! Maximum possible shift
      double precision bckgr       ! background

      double precision delta       ! Abstand Dekoration - zentrale Schicht
      double precision qa,qa2      ! Arguments
      double precision qhpt,qmod   ! fractions of qarg (mod Pi)
      double precision K           ! inner structure factor
      double precision K2,K2m1,Kn2 ! K^2, K^2-1, K^(n+2)
      double precision K3,K2p1     ! K^3, K^2+1
      double precision Keps
      double precision Ke1,Ke2
      double precision nn,nnn      ! n^2, n^3

      double precision A(31)       ! Coefficients of sin-oszillations
      double precision B(31)       ! arguments    of sin-oszillations

      integer          ier         ! Schultz distribution, error
      double precision az          ! Schultz distribution, Z+1
      double precision xpn         ! exponent, f_r=exp(xpn)

c     double precision DGAMMA      ! Gamma
c     double precision DLNGAMMA    ! lnGamma

      integer          KF
      double precision LGAMA

      double precision Oszi        ! function to calculate sum of sin-oszillations
      external         Oszi

      common /VesPar4/q,phi,no,nw,drho1,drho2,sig1,sig2,Rmain,zz,Rshift,
     x                bckgr

C--------------------------------------------


      VorF = phi * 4d-24*Pi*R0*R0 * 6d0/(((2d0*n+3d0)*n+1d0)*n)


      delta = 1.177410023d0 * (sig1+sig2)     ! form factor
      Fv    =  drho1*sig1*exp(-.5d0*q*q*sig1*sig1)
     +       + drho2*sig2*exp(-.5d0*q*q*sig2*sig2)*2d0*cos(delta*q)



      qa   = dabs(q*R0)                        ! structure factor
      qa2  = dabs(q*Rshift*(R0/Rmain))

      K = 3d0*(dsin(qa2)-qa2*dcos(qa2))/(qa2*qa2*qa2)    ! inner structure factor

      K2   = K*K
      K3   = K*K2
      K2m1 = K2 - 1d0
      K2p1 = K2 + 1d0
      Kn2  = K**(n+2d0)

      Keps = 1d-3
      Ke1  = Keps*420d0/36d0
     x       * ((((( 4d0*n+12d0)*n+13d0)*n+ 6d0)*n+ 1d0)*n*n)
     x       /((((((10d0*n+36d0)*n+21d0)*n-35d0)*n-35d0)*n*n+4d0)*n)
      Ke2  = Keps*0.5d0
     x       *  (((( 6d0*n+15d0)*n+10d0)*n*n-1d0)*n)
     x       /   ((( 6d0*n*n-10d0)*n*n+4d0)*n)

      qhpt= Pi*int(qa/Pi+0.5d0)               ! mod Pi
      qmod= qa - qhpt
      if (abs(qmod).lt.eps) then
        if (qmod.lt.0d0) then
          qmod = -eps
        else
          qmod =  eps
        endif
      endif
      qa  = qhpt + qmod

      nn  = n*n
      nnn = nn*n

      if (abs(K).le.Ke2) then

      Sv = ( n*(n+0.5d0)*(n+1d0)/6d0 + 0.25d0*
     x      (cos(2d0*qa*(n+1d0))*(-(n+1d0)/(sin(qa)**2d0)+(n+1d0)**2d0)
     x      +sin(2d0*qa*(n+1d0))*(-(n+1d0)**2d0+0.5d0/(sin(qa)**2d0) )
     x                          /tan(qa) ) ) / (qa*qa)

      else

      if (K.lt.1d0-Ke1) then

c2345678901234567890123456789012345678901234567890123456789012345678901234567xxxxxxxxx
      A(02)= -6d0*K2m1*K2p1*((K2+5d0)*K2+1d0)*nn
     &       +((((-6d0*K2-12d0)*K2+48d0)*K2+48d0)*K2+6d0)*n
     &       +(((( 3d0*K2+36d0)*K2+24d0)*K2-18d0)*K2-3d0)
      B(02)=qa*(2d0*n+1d0)

      A(03)=  (((( 3d0*K2-12d0)*K2-45d0)*K2-24d0)*K2-3d0)*nn
     &       +(((( 6d0*K2-12d0)*K2-72d0)*K2-48d0)*K2-6d0)*n
     &       +(((( 3d0*K2+18d0)*K2-24d0)*K2-36d0)*K2-3d0)
      B(03)=qa*(2d0*n-1d0)

      A(06)=  (((( 3d0*K2+24d0)*K2+45d0)*K2+12d0)*K2-3d0)*nn
     &       +6d0*K2*(3d0*K2+2d0)*n
     &       +3d0*K2*((4d0*K2-1d0)*K2-6d0)
      B(06)=qa*(2d0*n+3d0)

      A(07)= 18d0*K2*K2*K2p1*nn
     &       +6d0*K2*((6d0*K2+3d0)*K2-2d0)*n
     &       +3d0*K2*((6d0*K2+1d0)*K2-4d0)
      B(07)=qa*(2d0*n-3d0)

      A(10)=-18d0*K2*K2p1*nn
     &       -6d0*K2*K2*(2d0*K2+3d0)*n
     &       -3d0*K2*K2
      B(10)=qa*(2d0*n+5d0)

      A(11)=  3d0*K2*K2*nn
     &       +6d0*K2*K2*n
     &       +3d0*K2*K2
      B(11)=qa*(2d0*n-5d0)

      A(13)= -3d0*K2*K2*nn
      B(13)=qa*(2d0*n+7d0)


      A(26)=  6d0*K2m1*K2p1*((K2+4d0)*K2+1d0)*nnn
     &       +3d0*((K2-4d0)*K2-3d0)*((3d0*K2+8d0)*K2+1d0)*nn
     &       +(((( 3d0*K2-36d0)*K2-12d1)*K2-60d0)*K2-3d0)*n
     &      +42d0*K2*((K2-4d0)*K2+1d0)
      B(26)=qa

      A(28)= -2d0*K2m1*K2p1*(K2m1*K2+1d0)*nnn
     &       +((((-3d0*K2+ 9d0)*K2     )*K2+ 3d0)*K2+3d0)*nn
     &       +((((    -K2+ 7d0)*K2     )*K2+ 5d0)*K2+1d0)*n
     &       +6d0*K2*((-4d0*K2+11d0)*K2-4d0)
      B(28)=3d0*qa

      A(30)= -6d0*K2*K2m1*K2p1*nnn
     &       -3d0*K2*((K2-8d0)*K2-5d0)*nn
     &       +3d0*K2*((K2+8d0)*K2+3d0)*n
     &       +6d0*K2*( K2m1   *K2+1d0)
      B(30)=5d0*qa


      A(01)=  6d0*K3*((3d0*K2+10d0)*K2+5d0)*nn
     &       +6d0*K3*((3d0*K2+ 8d0)*K2+3d0)*n
     &      -12d0*K*K2m1*((K2+3d0)*K2+1d0)
      B(01)=2d0*qa*n

      A(05)=      K*(((-12d0*K2-12d0)*K2+12d0)*K2+ 6d0)*nn
     &       -6d0*K*(2d0*K2+3)*(2d0*K2m1*K2-1d0)*n
     &           +K*(((-12d0*K2p1   )*K2+36d0)*K2+12d0)
      B(05)=2d0*qa*(n-1d0)

      A(04)=      K* ((-30d0*K2-60d0)*K2-18d0)*nn
     &           +K* ((-42d0*K2-72d0)*K2-18d0)*n
     &           +K*(((-12d0*K2-36d0)*K2+12d0)*K2+12d0)
      B(04)=2d0*qa*(n+1d0)

      A(09)= -6d0*K3*(2*K2+1)*nn
     &       -6d0*K3*(4*K2+1)*n
     &      -12d0*K3*K2
      B(09)=2d0*qa*(n-2d0)

      A(08)= -6d0*K*(((      K2+ 2d0)*K2- 2d0)*K2- 2d0)*nn
     &       +6d0*K3*((K2+4d0)*K2+2d0)*n
     &      +12d0*K3
      B(08)=2d0*qa*(n+2d0)

      A(12)=  6d0*K3*(K2+2d0)*nn
     &       +6d0*K3*K2*n
      B(12)=2d0*qa*(n+3d0)


      A(27)= -6d0*K*K2m1*(2d0*K2+1d0)*(K2+2d0)*nnn
     &           +K*(((-12d0*K2+48d0)*K2+102d0)*K2+24d0)*nn
     &           +K* (( 66d0*K2+84d0)*K2+12d0)*n
     &      +24d0*K3*K2p1
      B(27)=2d0*qa

      A(29)=  6d0*K*K2m1*K2p1*K2p1*nnn
     &       +6d0*K*K2p1*((K2-5d0)*K2-2d0)*nn
     &       -6d0*K*K2p1*(5d0*K2+1d0)*n
     &      -12d0*K3*K2p1
      B(29)=4d0*qa

      A(31)=  2d0*K3*K2m1*nnn
     &       -6d0*K3*nn
     &       -2d0*K3*(K2+2d0)*n
      B(31)=6d0*qa


      A(14)=Kn2*K2m1*K*(-72d0             -12d0*n*(3d0*K2+4d0))
      B(14)=qa*n
      A(16)=Kn2*K2m1  *( 12d0*(3d0*K2-2d0)-12d0*n*(    K2+2d0))
      B(16)=qa*(n-1d0)
      A(15)=Kn2*K2m1  *(-12d0*(2d0*K2-3d0)+12d0*n*(4d0*K2+3d0))
      B(15)=qa*(n+1d0)
      A(18)=Kn2*K2m1*K*( 48d0             + 6d0*n*(4d0*K2+7d0))
      B(18)=qa*(n-2d0)
      A(17)=Kn2*K2m1*K*( 48d0             +12d0*n*(2d0*K2+1d0))
      B(17)=qa*(n+2d0)
      A(20)=Kn2*K2m1  *( -6d0*(4d0*K2-1d0)- 6d0*n*(2d0*K2-1d0))
      B(20)=qa*(n-3d0)
      A(19)=Kn2*K2m1  *(  6d0*(    K2-4d0)- 6d0*n*(7d0*K2+4d0))
      B(19)=qa*(n+3d0)
      A(22)=Kn2*K2m1*K*(-12d0             - 6d0*n*(    K2+2d0))
      B(22)=qa*(n-4d0)
      A(21)=Kn2*K2m1*K*(-12d0             - 6d0*n*(    K2-2d0))
      B(21)=qa*(n+4d0)
      A(24)=Kn2*K2m1*K2*6d0*(1d0+n)
      B(24)=qa*(n-5d0)
      A(23)=Kn2*K2m1  *(  6d0             + 6d0*n*(2d0*K2+1d0))
      B(23)=qa*(n+5d0)
      A(25)=Kn2*K2m1*K*(                  - 6d0*n)
      B(25)=qa*(n+6d0)


      Sv = Oszi(A,B)/(-48d0*sin(qa)**3d0*
     X                (K2p1-2d0*K*cos(qa))**4d0*qa*qa)

      else

      Sv = (-.50d0*cos(qa*(n+.5d0))*(n+.5d0)
     &      +.25d0*sin(qa*(n+.5d0))/tan(.5d0*qa))/(qa*sin(.5d0*qa))
      Sv = Sv*Sv

      endif

      endif

      if (zz.le.-.99d0) then

        f_r = 1d0

      else

        ier = -1                              ! Schultz distribution
        az = zz+1d0

        if (az.le.10d0) then
          KF = 1
          f_r = (az*R0/Rmain)**az/R0*dexp(-az*R0/Rmain)/LGAMA(KF,az)
        else
          KF = 0
          xpn = az*log(az)+zz*log(R0)-az*log(Rmain)-LGAMA(KF,az)
     X          -az*R0/Rmain
          xpn = min( 37d0,xpn)
          xpn = max(-37d0,xpn)
          f_r = dexp(xpn)
        endif

      endif

      inVesicle4 = (Vorf * 2d0*Pi * Fv*Fv * Sv  +  bckgr) * f_r

      end


c-------------------------------------------------------------------------------------------


      double precision function Oszi(A,B)

      implicit none

      double precision A(31)                  ! Amplitudes
      double precision B(31)                  ! Arguments

      double precision Pi
      parameter (Pi = 3.14159265358979323846264338328d0)

      integer          i,j,k                  ! loops
      integer          sgn                    ! sign
      integer          qitv                   ! interval
      double precision Sgn2                   ! sign
      double precision C(31)                  ! B(i) mod Pi/2
      integer          Case(31)               ! case sin/cos 0/1
      double precision Ca(31),Cs(31)          ! Cosine values (power of argument, sign)
      double precision Sa(31),Ss(31)          ! Sine   values (power of argument, sign)
      double precision T(20)                  ! Terms
      double precision fac                    ! faculty

c     double precision min                    ! minimum C
c     double precision KA,KB,KC,KCa,KCs,KSa,KSs ! copied values
c     integer          KCase

      do i=1,31

        sgn = 1
        if (B(i).lt.0d0) sgn = -1
        B(i) = dabs(B(i))

        qitv = int(B(i)*2d0/Pi+0.5d0)

        C(i) = B(i)-0.5d0*Pi*qitv             ! arg mod Pi/2

        Sgn2 = 1d0-2d0*(qitv/2-2*(qitv/4))    ! sign

        Case(i) = qitv-2*(qitv/2)             ! case sin/cos 0/1

        Ca(i) = 1d0
        Sa(i) = C(i)
        Cs(i) = Sgn2*sgn
        Ss(i) = Sgn2*sgn

      enddo


c      do i=1,30
c        k  = i
c        min= 1d2
c        do j=i,31
c          if(abs(C(j)).lt.min) then
c            k=j
c            min=abs(C(j))
c          endif
c        enddo
c        KA = A(i)
c        KB = B(i)
c        KC = C(i)
c        KCa= Ca(i)
c        KCs= Cs(i)
c        KSa= Sa(i)
c        KSs= Ss(i)
c        KCase = Case(i)
c        A(i) = A(k)
c        B(i) = B(k)
c        C(i) = C(k)
c        Ca(i)= Ca(k)
c        Cs(i)= Cs(k)
c        Sa(i)= Sa(k)
c        Ss(i)= Ss(k)
c        Case(i) = Case(k)
c        A(k) = KA
c        B(k) = KB
c        C(k) = KC
c        Ca(k)= KCa
c        Cs(k)= KCs
c        Sa(k)= KSa
c        Ss(k)= KSs
c        Case(k) = KCase
c      enddo

      fac = 1d0
      j = 1

      do k=0,9

        T(j) = 0d0
        do i=1,31
          if(Case(i).eq.1) then
            T(j) = T(j) + Cs(i)*Ca(i)*A(i)
            Ca(i)= Ca(i)*(C(i)*C(i))
            Cs(i)=-Cs(i)
          endif
        enddo
        T(j)=T(j)/fac
        fac =fac*j
        j=j+1

        T(j) = 0d0
        do i=1,31
          if(Case(i).eq.0) then
            T(j) = T(j) + Ss(i)*Sa(i)*A(i)
            Sa(i)= Sa(i)*(C(i)*C(i))
            Ss(i)=-Ss(i)
          endif
        enddo
        T(j)=T(j)/fac
        fac =fac*j
        j=j+1

      enddo

      Oszi = 0d0

      do j=1,20
        Oszi = Oszi + T(j)
      enddo

      end



        double precision function LGAMA(KF,X)
C
C       ==================================================
C       Purpose: Compute gamma function â(x) or ln[â(x)]
C       Input:   x  --- Argument of â(x) ( x > 0 )
C                KF --- Function code
C                       KF=1 for â(x); KF=0 for ln[â(x)]
C       Output:  GL --- â(x) or ln[â(x)]
C       ==================================================
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION A(10)
        DATA A/8.333333333333333D-02,-2.777777777777778D-03,
     &         7.936507936507937D-04,-5.952380952380952D-04,
     &         8.417508417508418D-04,-1.917526917526918D-03,
     &         6.410256410256410D-03,-2.955065359477124D-02,
     &         1.796443723688307D-01,-1.39243221690590D+00/
        X0=X
        IF (X.EQ.1.0.OR.X.EQ.2.0) THEN
           GL=0.0D0
           GO TO 20
        ELSE IF (X.LE.7.0) THEN
           N=INT(7-X)
           X0=X+N
        ENDIF
        X2=1.0D0/(X0*X0)
        XP=6.283185307179586477D0
        GL0=A(10)
        DO 10 K=9,1,-1
10         GL0=GL0*X2+A(K)
        GL=GL0/X0+0.5D0*DLOG(XP)+(X0-.5D0)*DLOG(X0)-X0
        IF (X.LE.7.0) THEN
           DO 15 K=1,N
              GL=GL-DLOG(X0-1.0D0)
15            X0=X0-1.0D0
        ENDIF
20      IF (KF.EQ.1) GL=DEXP(GL)
        LGAMA = GL
        RETURN
        END

