C     =================================================================
C     File: film.f
C     =================================================================

C     =================================================================
C     Module: Subroutines that define the problem
C     =================================================================

C     Last update of any of the component of this module: 
 
C     April 30, 2008.

C     =================================================================

C     Film problem
C     -------------
C
C     This problem consists of finding eta and kappa in R^np such that 
C     the distance between the trasmitance observed t^{obs} and the
C     transmitance calculated T(lambda_i,s,d,eta_i,kappa_i) (lambda_i, 
C     s and d given) is as small as possible:
C
C     Minimize sum_{i = 1}^{np} [T^{obs}(lambda_i) - 
C                                T(lambda_i,s,d,eta_i,kappa_i)]^2
C 
C     subject to   eta_{i+1} <= eta_i,   i = 1, ..., np-1,
C                kappa_{i+1} <= kappa_i, i = 1, ..., np-1,
C                      eta_i <= eta_{i-1} + (eta_{i+1} - eta_{i-1})*
C            (lambda_i - lambda_{i-1})/(lambda_{i+1} - lambda_{i-1}), 
C                                                     i = 2, ..., np-1,
C
C                    kappa_i >= kappa_{i-1} + (kappa_{i+1} - kappa_{i-1})*
C            (lambda_i - lambda_{i-1})/(lambda_{i+1} - lambda_{i-1}), 
C                                        lambda_{i+1} <= lambda_{infl},
C
C                    kappa_i <= kappa_{i-1} + (kappa_{i+1} - kappa_{i-1})*
C            (lambda_i - lambda_{i-1})/(lambda_{i+1} - lambda_{i-1}), 
C                                        lambda_{i-1} >= lambda_{infl},
C
C                      eta_i >= 1, i = 1, ..., np
C                    kappa_i >= 0, i = 1, ..., np
C
C     where eta_i = eta(lambda_i) and kappa_i = kappa(lambda_i).

C     ******************************************************************
C     ******************************************************************

      subroutine inip(n,x,l,u,ml,mleq,lda,A,b,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer n,ml,mleq,lda,inform

C     ARRAY ARGUMENTS
      double precision l(*),u(*),x(*),A(lda,*),b(*)

C     This subroutine must set some problem data. For achieving this 
C     objective YOU MUST MODIFY it according to your problem. See below 
C     where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     This subroutine has no input parameters.
C
C     On Return
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              initial point,
C
C     l        double precision l(n),
C              lower bounds on x,
C
C     u        double precision u(n),
C              upper bounds on x,
C
C     ml       integer,
C              total number of linear constraints (excluding the bounds),
C
C     mleq     integer,
C              number of linear equality constraints (excluding the bounds),
C
C     lda      integer,
C              leading dimension of matrix A,
C
C     A        double precision A(ml,n)
C              coeficient matrix of linear constraints, such that Ax = b, 
C              for rows 1 to mleq of A and Ax >= b, for rows mleq+1 to ml,
C
C     b        double precision b(ml)
C              right-hand side of linear constraints,

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR PROBLEM
C     DATA:
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued,de1,de2

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i,j
      double precision tmp,dkick

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

      inform = 0

C     Generate problem

      call genprob(dkick)
      d          = trued
      indlaminfl = 1

C     Read problem data

      open(12, file="prob.dat")

      do i = 1, np
         read(12,*) tmp
         tobs(i) = tmp
      end do

      close(12)

C     Set initial point

      de1 = 4.0d0
      de2 = 2.0d0
      call geninp(np,lambmin,lambmax,de1,de2,x)

C     Set number of variables

      n = 2 * np

C     First np components of x represent eta and last np components of 
C     x respresent kappa

C     Lower and upper bounds

C     eta_i >= 1

      do i = 1,np
          l(i) = 1.0d+00
          u(i) = 1.0d+20
      end do

C     kappa_i >= 0

      do i = np+1,n
          l(i) = 0.0d+00
          u(i) = 1.0d+20
      end do

C     Number of linear constraints (equalities plus inequalities)

      ml   = 4 * np - 6
      mleq = 0
 
C     Set linear constraints

      do j = 1,n
         do i = 1,ml
            A(i,j) = 0.0d0
         end do
      end do

      do i = 1,ml
         b(i) = 0.0d0
      end do

      j = 0

C     eta_{i+1} <= eta_i (eta_1 to eta_np stored in x(1) to x(np))

      do i = 1,np-1
         j = j + 1
         A(j,i)   =  1.0d0
         A(j,i+1) = -1.0d0
      end do

C     kappa_{i+1} <= kappa_i (kappa_1 to kappa_np stored in x(np+1) to x(2*np))

      do i = 1,np-1
         j = j + 1
         A(j,np+i)   =  1.0d0
         A(j,np+i+1) = -1.0d0
      end do

C     eta_i <= eta_{i-1} + 
C     (eta_{i+1} - eta{i-1})/(lamb_{i+1} - lamb_{i-1}) * 
C     (lamb_i - lamb_{i-1}), with
C     lamb_i = lambmin + (i-1)*(lambmax-lambmin)/(np-1)

      do i = 2,np-1
         j = j + 1
         A(j,i-1) =  0.5d0
         A(j,i)   = -1.0d0 
         A(j,i+1) =  0.5d0
      end do

C     kappa_i >= kappa_{i-1} + 
C     (kappa_{i+1} - kappa{i-1})/(lamb_{i+1} - lamb_{i-1}) * 
C     (lamb_i - lamb_{i-1}), for i = 2 to indlaminfl, with
C     lamb_i = lambmin + (i-1)*(lambmax-lambmin)/(np-1)

      do i = 2,indlaminfl
         j = j + 1
         A(j,np+i-1) = -0.5d0
         A(j,np+i)   =  1.0d0 
         A(j,np+i+1) = -0.5d0
      end do

C     kappa_i <= kappa_{i-1} + 
C     (kappa_{i+1} - kappa{i-1})/(lamb_{i+1} - lamb_{i-1}) * 
C     (lamb_i - lamb_{i-1}), for i = indlaminfl+1 to np-1, with
C     lamb_i = lambmin + (i-1)*(lambmax-lambmin)/(np-1)

      do i = indlaminfl+1,np-1
         j = j + 1
         A(j,np+i-1) =  0.5d0
         A(j,np+i)   = -1.0d0 
         A(j,np+i+1) =  0.5d0
      end do

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE INIP.
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalf(n,x,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     This subroutine must compute the objective function. For achieving 
C     this objective YOU MUST MODIFY it according to your problem. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     f        double precision,
C              objective function value at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the objective 
C              function. (For example, trying to compute the square root 
C              of a negative number, dividing by zero or a very small 
C              number, etc.) If everything was o.k. you must set it 
C              equal to zero.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR OBJECTIVE
C     FUNCTION:
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision tr(npmax)

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

      flag = 0

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr

      call compt(np,lambmin,lambmax,s,d,x,tr)

      f = 0.0d0
      do i = 1,np
         f = f + (tobs(i)-tr(i))**2
      end do

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALF.
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalg(n,x,g,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n

C     ARRAY ARGUMENTS
      double precision g(n),x(n)

C     This subroutine must compute the gradient vector of the objective 
C     function. For achieving these objective YOU MUST MODIFY it in the 
C     way specified below. However, if you decide to use numerical 
C     derivatives (we dont encourage this option at all!) you dont need
C     to modify evalg.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     g        double precision g(n),
C              gradient vector of the objective function evaluated at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of any component 
C              of the gradient vector. (For example, trying to compute 
C              the square root of a negative number, dividing by zero or 
C              a very small number, etc.) If everything was o.k. you 
C              must set it equal to zero.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENT
C     VECTOR OF YOUR OBJECTIVE FUNCTION: 
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision tr(npmax),dtr(2*npmax)

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

      flag = 0

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr and
C            dT[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C            dT[lamb_i,s_i,d,eta_i,kappa_i] = dtr, where
C            dtr(i) first derivative of T with respect to eta_i, 
C                    for i = 1 to np and
C            dtr(i) first derivative of T with respect to kappa_i, 
C                    for i = np+1 to 2*np

      call compdt(np,lambmin,lambmax,s,d,x,tr,dtr)

      do i = 1,np
         g(i) = - 2.0d0 * (tobs(i)-tr(i)) * dtr(i)
      end do

      do i = np+1,n
         g(i) = - 2.0d0 * (tobs(i-np)-tr(i-np)) * dtr(i)
      end do

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALG. 
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalh(n,x,ldh,H,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,ldh

C     ARRAY ARGUMENTS
      double precision x(n),H(ldh,n)

C     This subroutine might compute the Hessian matrix of the objective 
C     function. For achieving this objective YOU MAY MODIFY it according 
C     to your problem. To modify this subroutine IS NOT MANDATORY. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     ldh      integer,
C              leading dimension of matrix H,
C
C     On Return
C
C     H        double precision H(ldh,n),
C              Hessian matrix of the objective function. Just the lower 
C              triangular part of Hessian matrix must be computed,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the Hessian
C              matrix of the objective funtion. (For example, trying 
C              to compute the square root of a negative number, 
C              dividing by zero or a very small number, etc.) If 
C              everything was o.k. you must set it equal to zero.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE HESSIAN
C     MATRIX OF YOUR OBJECTIVE FUNCTION (DENSE FORMAT): 
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision tr(npmax),dtr(2*npmax)

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr;
C            dT[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C            dT[lamb_i,s_i,d,eta_i,kappa_i] = dtr, where
C            dtr(i) first derivative of T with respect to eta_i, 
C                   for i = 1 to np and
C            dtr(i) first derivative of T with respect to kappa_i, 
C                   for i = np+1 to 2*np
C    and
C           d2T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C           d2T[lamb_i,s_i,d,eta_i,kappa_i] = Htr, where
C           Htr(i,j) second derivative of T with respect to eta_i 
C                   and eta_j, for i = 1 to np, j = 1 to np
C           Htr(i,j) second derivative of T with respect to kappa_i 
C                   and kappa_j for i = np+1 to 2*np, j = np+1 to 2*np
C           Htr(i,j) second derivative of T with respect to eta_i 
C                   and kappa_j for i = 1 to np, j = np+1 to 2*np

      flag = 0

      call compd2t(np,lambmin,lambmax,s,d,x,tr,dtr,ldh,H)

      do i = 1, np
         H(i,i) = 2.0d0*dtr(i)**2 - 2.0d0*(tobs(i)-tr(i))*H(i,i)
      end do

      do i = np+1, n
         H(i,i) = 2.0d0*dtr(i)**2 - 2.0d0*(tobs(i-np)-tr(i-np))*H(i,i)
      end do

      do i = 1, np
         H(i+np,i) = - 2.0d0 * (tobs(i)-tr(i)) * H(i+np,i)
      end do

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALH. 
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalsh(n,x,hlin,hcol,hval,hnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*)
      double precision hval(*),x(n)

C     This subroutine might compute the Hessian matrix of the objective 
C     function. For achieving this objective YOU MAY MODIFY it according 
C     to your problem. To modify this subroutine IS NOT MANDATORY. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     hnnz     integer,
C              number of perhaps-non-null elements of the computed 
C              Hessian,
C
C     hlin     integer hlin(hnnz),
C              see below,
C
C     hcol     integer hcol(hnnz),
C              see below,
C
C     hval     double precision hval(hnnz),
C              the non-null value of the (hlin(k),hcol(k)) position 
C              of the Hessian matrix of the objective function must 
C              be saved at hval(k). Just the lower triangular part of
C              Hessian matrix must be computed,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the Hessian
C              matrix of the objective funtion. (For example, trying 
C              to compute the square root of a negative number, 
C              dividing by zero or a very small number, etc.) If 
C              everything was o.k. you must set it equal to zero.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE HESSIAN
C     MATRIX OF YOUR OBJECTIVE FUNCTION (SPARSE FORMAT): 
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i,j,lin,col
      double precision val

C     LOCAL ARRAYS
      double precision tr(npmax),dtr(2*npmax)

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr;
C            dT[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C            dT[lamb_i,s_i,d,eta_i,kappa_i] = dtr, where
C            dtr(i) first derivative of T with respect to eta_i, 
C                   for i = 1 to np and
C            dtr(i) first derivative of T with respect to kappa_i, 
C                   for i = np+1 to 2*np
C    and
C           d2T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C           d2T[lamb_i,s_i,d,eta_i,kappa_i] = Htr, where
C           Htr(i,j) second derivative of T with respect to eta_i 
C                   and eta_j, for i = 1 to np, j = 1 to np
C           Htr(i,j) second derivative of T with respect to kappa_i 
C                   and kappa_j for i = np+1 to 2*np, j = np+1 to 2*np
C           Htr(i,j) second derivative of T with respect to eta_i 
C                   and kappa_j for i = 1 to np, j = np+1 to 2*np

      flag = 0

      call compsd2t(np,lambmin,lambmax,s,d,x,tr,dtr,hnnz,hlin,hcol,hval)

      do i = 1, hnnz
         lin = hlin(i)
         col = hcol(i)
         val = hval(i)

         if (lin .eq. col) then 
            if (lin .gt. np) then
               j = lin - np
            else
               j = lin
            end if
            hval(i) = 2.0d0*dtr(j)**2 - 2.0d0*(tobs(j)-tr(j))*val
         else
            j = min(lin,col)
            hval(i) = - 2.0d0 * (tobs(j)-tr(j)) * val
         end if
      end do

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALSH. 
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhp(n,x,p,hp,goth,flag)

      implicit none

C     SCALAR ARGUMENTS
      logical goth
      integer flag,n

C     ARRAY ARGUMENTS
      double precision hp(n),p(n),x(n)

C     This subroutine might compute the product of the Hessian of the
C     Lagrangian times vector p (just the Hessian of the objective 
C     function in the unconstrained or bound-constrained case). 
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     m        integer,
C              number of constraints,
C
C     lambda   double precision lambda(m),
C              vector of Lagrange multipliers,
C
C     p        double precision p(n),
C              vector of the matrix-vector product,
C
C     goth     logical,
C              can be used to indicate if the Hessian matrices were
C              computed at the current point. It is set to .false.
C              by the optimization method every time the current
C              point is modified. Sugestion: if its value is .false. 
C              then compute the Hessians, save them in a common 
C              structure and set goth to .true.. Otherwise, just use 
C              the Hessians saved in the common block structure,
C
C     On Return
C
C     hp       double precision hp(n),
C              Hessian-vector product,
C
C     goth     logical,
C              see above,
C              
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the 
C              Hessian-vector product. (For example, trying to compute 
C              the square root of a negative number, dividing by zero 
C              or a very small number, etc.) If everything was o.k. you 
C              must set it equal to zero.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE PRODUCT OF
C     HESSIAN MATRIX OF YOUR OBJECTIVE FUNCTION AND A VECTOR P:
C     ******************************************************************

      flag = -1

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALHP. 
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************

      subroutine endp(n,x,l,u,ml,mleq,lda,A,b)

      implicit none

C     SCALAR ARGUMENTS
      integer n,ml,mleq,lda

C     ARRAY ARGUMENTS
      double precision l(n),u(n),x(n),A(lda,n),b(ml)

C     This subroutine can be used to do some extra job after the solver
C     has found the solution,like some extra statistics, or to save the
C     solution in some special format or to draw some graphical
C     representation of the solution. If the information given by the
C     solver is enough for you then leave the body of this subroutine
C     empty.
C     
C     Parameters of the subroutine:
C
C     The parameters of this subroutine are the same parameters of
C     subroutine inip. But in this subroutine there are not output
C     parameter. All the parameters are input parameters.

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SOME EXTRA JOB YOU
C     WANT AFTER THE METHOD IS EXECUTED:
C     ******************************************************************

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i
      double precision ei, ai

C     LOCAL ARRAYS
      double precision tr(npmax),lamb(npmax)

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

      call compt(np,lambmin,lambmax,s,d,x,tr)

      do i = 1, np
         lamb(i) = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
      end do

C     Write lambda X transmition (true and retrieved)

      open(21, file='transmition.out')

      do i = 1, np
         write(21,100) lamb(i), tobs(i), tr(i)
      end do

      close(21)

C     Write lambda X eta (true and retrieved)

      open(22, file='refraction.out')

      do i = 1, np
         write(22,100) lamb(i), teta(i), x(i)
      end do

      close(22)

C     Write lambda X eta (true and retrieved)

      open(23, file='attenuation.out')

      do i = 1, np
         write(23,100) lamb(i), tkappa(i), x(np+i)
      end do

      close(23)

C     Write energy X eta (true and retrieved)

      open(24, file='absorption.out')

      do i = 1, np
         ei = 1.24d3/lamb(i)
         ai = (1.0d7 * 4.0d0 * pi * x(np+i))/lamb(i) 
         write(24,100) ei, log10(1.0d7*talpha(i)), log10(ai)
      end do

      close(24)

 100  format(1X,F20.6,1X,F20.6,1X,F20.6)


C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE ENDP. 
C     ******************************************************************

      end

C     ******************************************************************
C     FROM THIS POINT, ALL THE SUBROUTINE ARE AUXILIARY. ALL EXTRA 
C     SUBROUTINE YOU NEED MUST BE PUT HERE AND THESE ONES MIGHT BE 
C     DELETED WHEN YOU SET YOUR OWN PROBLEM.
C     ******************************************************************

C     ******************************************************************
C     ******************************************************************

      subroutine compt(np,lambmin,lambmax,s,d,x,tr)

      implicit none

C     PARAMETERS
      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     SCALAR ARGUMENTS
      integer np
      double precision lambmin,lambmax,d

C     ARRAY ARGUMENTS
      double precision s(np),x(2*np),tr(np)

C     LOCAL SCALARS
      integer i
      double precision eta,eta2,kappa,kappa2,s2,lambi,phi,xx,xx2,AA,B1,
     +     B2,BB,C1,C2,C3,C4,C5,C6,CC,D1,D2,DD

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr

      do i = 1, np
         eta    = x(i)
         eta2   = eta**2
         kappa  = x(i+np)
         kappa2 = kappa**2
         s2     = s(i)**2
         
         lambi  = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
         phi    = (4.0d0 * pi * d * eta)/lambi
         xx     = exp((-4.0d0 * pi * kappa * d)/lambi)
         xx2    = xx**2

         AA     = 16.0d0 * s(i) * (eta2 + kappa2)

         B1     = (eta + 1.0d0)**2 + kappa2
         B2     = (eta + 1.0d0) * (eta + s2) + kappa2
         BB     = B1 * B2

         C1     = (eta2 + kappa2 - 1.0d0)
         C2     = (eta2 + kappa2 - s2)
         C3     = 2.0d0 * kappa2 * (s2 + 1.0d0)
         C4     = 2.0d0 * cos(phi)
         C5     = (s2 + 1.0d0) * C1
         C6     = 2.0d0 * sin(phi) * kappa
         CC     = (C1 * C2 - C3) * C4 - (2.0d0 * C2 + C5) * C6

         D1     = (eta - 1.0d0)**2 + kappa2
         D2     = (eta - 1.0d0) * (eta - s2) + kappa2
         DD     = D1 * D2

         tr(i)   = (AA * xx) / (BB - CC * xx + DD * xx2)
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compdt(np,lambmin,lambmax,s,d,x,tr,dtr)

      implicit none

C     PARAMETERS
      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     SCALAR ARGUMENTS
      integer np
      double precision lambmin,lambmax,d

C     ARRAY ARGUMENTS
      double precision s(np),x(2*np),tr(np),dtr(2*np)

C     LOCAL SCALARS
      integer i
      double precision eta,eta2,kappa,kappa2,s2,lambi,phi,xx,xx2,AA,B1,
     +     B2,BB,C1,C2,C3,C4,C5,C6,CC,D1,D2,DD,Ax,beta,dAA,dB1,dB2,dBB,
     +     dC1,dC2,dC3,dC4,dC5,dC6,dCC,dD1,dD2,dDD,dAx,dbeta,dx

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr
 
      do i = 1, np
         eta    = x(i)
         eta2   = eta**2
         kappa  = x(i+np)
         kappa2 = kappa**2
         s2     = s(i)**2
         
         lambi  = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
         phi    = (4.0d0 * pi * d * eta)/lambi
         xx     = exp((-4.0d0 * pi * kappa * d)/lambi)
         xx2    = xx**2

         AA     = 16.0d0 * s(i) * (eta2 + kappa2)

         B1     = (eta + 1.0d0)**2 + kappa2
         B2     = (eta + 1.0d0) * (eta + s2) + kappa2
         BB     = B1 * B2

         C1     = (eta2 + kappa2 - 1.0d0)
         C2     = (eta2 + kappa2 - s2)
         C3     = 2.0d0 * kappa2 * (s2 + 1.0d0)
         C4     = 2.0d0 * cos(phi)
         C5     = (s2 + 1.0d0) * C1
         C6     = 2.0d0 * sin(phi) * kappa
         CC     = (C1 * C2 - C3) * C4 - (2.0d0 * C2 + C5) * C6

         D1     = (eta - 1.0d0)**2 + kappa2
         D2     = (eta - 1.0d0) * (eta - s2) + kappa2
         DD     = D1 * D2

         tr(i)   = (AA * xx) / (BB - CC * xx + DD * xx2)

C        Derivatives with respect to eta

         Ax     = AA * xx
         beta   = (BB - CC * xx + DD * xx2)

         dAA    = 32.0d0 * s(i) * eta

         dB1    = 2.0d0 * (eta + 1.0d0)
         dB2    = 2.0d0 * eta + s2 + 1.0d0
         dBB    = B1 * dB2 + dB1 * B2

         dC1    = 2.0d0 * eta
         dC2    = 2.0d0 * eta
         dC3    = 0.0d0
         dC4    = - (8.0d0 * sin(phi) * pi * d )/lambi
         dC5    = 4.0d0 * eta
         dC5    = 2.0d0 * (s2 + 1.0d0) * eta
         dC6    = (8.0d0 * kappa * cos(phi) * d * pi)/lambi
         dCC    = (C1 * C2 - C3) * dC4 + (C1 * dC2 + dC1 * C2) * C4 
     +          - (2.0d0 * C2 + C5)*dC6 - (2.0d0 * dC2 + dC5) * C6

         dD1    = 2.0d0 * (eta - 1.0d0)
         dD2    = 2.0d0 * eta - s2 - 1.0d0
         dDD    = D1 * dD2 + dD1 * D2

         dx     = 0.0d0

         dAx    = dAA * xx
         dbeta  = dBB - dCC * xx + dDD * xx2

         dtr(i) = (dAx * beta - Ax * dbeta) / beta**2

C        Derivatives with respect to kappa

         dAA    = 32.0d0 * s(i) * kappa

         dB1    = 2.0d0 * kappa
         dB2    = 2.0d0 * kappa
         dBB    = B1 * dB2 + dB1 * B2

         dC1    = 2.0d0 * kappa
         dC2    = 2.0d0 * kappa
         dC3    = 4.0d0 * kappa * (s2 + 1.0d0)
         dC4    = 0.0d0
         dC5    = 4.0d0 * kappa
         dC5    = 2.0d0 * (s2 + 1.0d0) * kappa
         dC6    = 2.0d0 * sin(phi)
         dCC    = (C1 * dC2 + dC1 * C2 - dC3) * C4 - (2.0d0*C2 + C5)*dC6 
     +          - (2.0d0 * dC2 + dC5) * C6

         dD1    = 2.0d0 * kappa
         dD2    = 2.0d0 * kappa
         dDD    = D1 * dD2 + dD1 * D2

         dx     = (-4.0d0 * pi * d * xx)/lambi

         dAx    = dAA * xx + AA * dx
         dbeta  = dBB - (dCC*xx + CC*dx) + (dDD*xx2 + DD * 2.0d0*xx*dx)

         dtr(i+np) = (dAx * beta - Ax * dbeta) / beta**2
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compd2t(np,lambmin,lambmax,s,d,x,tr,dtr,ldhtr,Htr)

      implicit none

C     PARAMETERS
      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     SCALAR ARGUMENTS
      integer np,ldhtr
      double precision lambmin,lambmax,d

C     ARRAY ARGUMENTS
      double precision s(np),x(2*np),tr(np),dtr(2*np),Htr(ldhtr,2*np)

C     LOCAL SCALARS
      integer i,j,n
      double precision eta,eta2,kappa,kappa2,s2,lambi,phi,xx,xx2,tmp,AA,
     +     B1,B2,BB,C1,C2,C3,C4,C5,C6,CC,D1,D2,DD,Ax,beta,dAAde,dB1de,
     +     dB2de,dBBde,dC1de,dC2de,dC4de,dC5de,dC6de,dCCde,dD1de,dD2de,
     +     dDDde,dAxde,dbetade,dAAdk,dB1dk,dB2dk,dBBdk,dC1dk,dC2dk,
     +     dC3dk,dC5dk,dC6dk,dCCdk,dD1dk,dD2dk,dDDdk,dAxdk,dbetadk,dxdk,
     +     dAAdede,dBBdede,dC4dede,dC6dede,dCCdede,dDDdede,dAxdede,
     +     dbetadede,dBBdedk,dC6dedk,dCCdedk,dDDdedk,dAxdedk,dbetadedk,
     +     dAAdkdk,dBBdkdk,dC3dkdk,dCCdkdk,dDDdkdk,dAxdkdk,dbetadkdk

      n = 2 * np

      do j = 1,n
         do i = j,n
            Htr(i,j) = 0.0d0
         end do
      end do

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr
 
      do i = 1, np
         eta    = x(i)
         eta2   = eta**2
         kappa  = x(i+np)
         kappa2 = kappa**2
         s2     = s(i)**2
         
         lambi  = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
         phi    = (4.0d0 * pi * d * eta)/lambi
         xx     = exp((-4.0d0 * pi * kappa * d)/lambi)
         xx2    = xx**2
         
         AA     = 16.0d0 * s(i) * (eta2 + kappa2)
         
         B1     = (eta + 1.0d0)**2 + kappa2
         B2     = (eta + 1.0d0) * (eta + s2) + kappa2
         BB     = B1 * B2
         
         C1     = (eta2 + kappa2 - 1.0d0)
         C2     = (eta2 + kappa2 - s2)
         C3     = 2.0d0 * kappa2 * (s2 + 1.0d0)
         C4     = 2.0d0 * cos(phi)
         C5     = (s2 + 1.0d0) * C1
         C6     = 2.0d0 * sin(phi) * kappa
         CC     = (C1 * C2 - C3) * C4 - (2.0d0 * C2 + C5) * C6
         
         D1     = (eta - 1.0d0)**2 + kappa2
         D2     = (eta - 1.0d0) * (eta - s2) + kappa2
         DD     = D1 * D2
         
         tr(i)   = (AA * xx) / (BB - CC * xx + DD * xx2)
         
         Ax      = AA * xx
         beta    = (BB - CC * xx + DD * xx2)
         
C        First derivatives with respect to eta
            
         dAAde   = 32.0d0 * s(i) * eta
         
         dB1de   = 2.0d0 * (eta + 1.0d0)
         dB2de   = 2.0d0 * eta + s2 + 1.0d0
         dBBde   = B1 * dB2de + dB1de * B2
         
         dC1de   = 2.0d0 * eta
         dC2de   = 2.0d0 * eta
         dC4de   = - (8.0d0 * sin(phi) * pi * d )/lambi
         dC5de   = 4.0d0 * eta
         dC6de   = (8.0d0 * kappa * cos(phi) * d * pi)/lambi
         dCCde   = (C1 * C2 - C3) * dC4de + (C1*dC2de + dC1de*C2) * C4 
     +           - (2.0d0*C2 + C5) * dC6de - (2.0d0*dC2de + dC5de) * C6
         
         dD1de   = 2.0d0 * (eta - 1.0d0)
         dD2de   = 2.0d0 * eta - s2 - 1.0d0
         dDDde   = D1 * dD2de + dD1de * D2
         
         dAxde   = dAAde * xx
         dbetade = dBBde - dCCde * xx + dDDde * xx2
         
         dtr(i) = (dAxde * beta - Ax * dbetade) / beta**2
         
C        First derivatives with respect to kappa

         dAAdk   = 32.0d0 * s(i) * kappa
         
         dB1dk   = 2.0d0 * kappa
         dB2dk   = 2.0d0 * kappa
         dBBdk   = B1 * dB2dk + dB1dk * B2
         
         dC1dk   = 2.0d0 * kappa
         dC2dk   = 2.0d0 * kappa
         dC3dk   = 4.0d0 * kappa * (s2 + 1.0d0)
         dC5dk   = 4.0d0 * kappa
         dC6dk   = 2.0d0 * sin(phi)
         dCCdk   = (C1 * dC2dk + dC1dk * C2 - dC3dk) * C4 
     +           - (2.0d0*C2 + C5) * dC6dk - (2.0d0*dC2dk + dC5dk) * C6
         
         dD1dk   = 2.0d0 * kappa
         dD2dk   = 2.0d0 * kappa
         dDDdk   = D1 * dD2dk + dD1dk * D2
         
         dxdk    = (-4.0d0 * pi * d * xx)/lambi
         
         dAxdk   = dAAdk * xx + AA * dxdk
         dbetadk = dBBdk - (dCCdk*xx + CC*dxdk) 
     +           + (dDDdk*xx2 + DD*2.0d0*xx*dxdk)
         
         dtr(i+np) = (dAxdk * beta - Ax * dbetadk) / beta**2
         
C        Second derivatives with respect to eta and eta

         dAAdede   = 32.0d0 * s(i)
         
         dBBdede   = 2.0d0 * (dB1de * dB2de + B1 + B2)
         
         dC4dede   = - (32.0d0 * cos(phi) * pi**2 * d**2 )/lambi**2
         dC6dede   = -(32.0d0*kappa*sin(phi)*d**2*pi**2)/lambi**2
         dCCdede   = (dC1de * C2 + C1 * dC2de) * dC4de 
     +             + (C1*C2 - C3) * dC4dede 
     +             + 2.0d0*C4*(dC1de * dC2de + C1 + C2) 
     +             + dC4de * (C1 * dC2de + dC1de * C2) 
     +             - 2.0d0 * dC6de * (2.0d0*dC2de+dC5de) 
     +             - (2.0d0*C2 + C5) * dC6dede - 8.0d0 * C6
         
         dDDdede   = 2.0d0 * (dD1de * dD2de + D1 + D2)
         
         dAxdede   = dAAdede * xx 
         dbetadede = dBBdede - dCCdede*xx + dDDdede*xx2
            
         tmp = (dAxdede*beta - Ax*dbetadede) * beta**2 
     +       - 2.0d0*beta*dbetade * (dAxde*beta - Ax*dbetade)
         Htr(i,i) = tmp/(beta**4)
         
C        Second derivatives with respect to eta and kappa

         dBBdedk   = dB1dk * dB2de + dB1de * dB2dk
         
         dC6dedk   = (8.0d0 * cos(phi) * d * pi)/lambi
         dCCdedk   = (dC1dk * C2 + C1 * dC2dk - dC3dk) * dC4de 
     +             + (dC1dk * dC2de + dC1de * dC2dk) * C4 
     +             - (2.0d0 * dC2dk + dC5dk) * dC6de 
     +             - (2.0d0 * C2 + C5) * dC6dedk 
     +             - (2.0d0 * dC2de + dC5de) * dC6dk
         
         dDDdedk   = dD1dk * dD2de + dD1de * dD2dk
         
         dAxdedk   = dAAde * dxdk
         dbetadedk = dBBdedk - dCCdedk*xx - dCCde*dxdk + dDDdedk*xx2 
     +             + 2.0d0*dDDde*xx*dxdk
         
         tmp = (dAxdedk*beta + dAxde*dbetadk - dAxdk*dbetade 
     +        - Ax*dbetadedk) * beta**2 - 2.0d0*beta*dbetadk 
     +        * (dAxde*beta - Ax*dbetade)
         
         Htr(i+np,i) = tmp/(beta**4)
         
C        Second derivatives with respect to kappa and kappa
         
         dAAdkdk   = 32.0d0 * s(i)
         
         dBBdkdk   = 2.0d0 * (dB1dk*dB2dk + B1 + B2)
         
         dC3dkdk   = 4.0d0 * (s2 + 1.0d0)
         dCCdkdk   = 2.0d0*C4*(dC1dk * dC2dk + C1 + C2) - dC3dkdk * C4 
     +             - 2.0d0 * dC6dk * (2.0d0*dC2dk + dC5dk) - 8.0d0 * C6
         
         dDDdkdk   = 2.0d0 * (dD1dk*dD2dk + D1 + D2)
         
         dAxdkdk   = dAAdkdk * xx + 2.0d0 * dAAdk * dxdk
         dbetadkdk = dBBdkdk - dCCdkdk*xx - 2.0d0*dCCdk*dxdk+dDDdkdk*xx2 
     +             + 2.0d0*dDDdk*xx*dxdk + 2.0d0*dxdk*(dDDdk*xx+DD*dxdk)
         
         tmp = (dAxdkdk*beta - Ax*dbetadkdk) * beta**2 
     +       - 2.0d0*beta*dbetadk * (dAxdk*beta - Ax*dbetadk)
         Htr(i+np,i+np) = tmp/(beta**4)
      end do
      
      end

C     ******************************************************************
C     ******************************************************************

      subroutine compsd2t(np,lambmin,lambmax,s,d,x,tr,dtr,hnnz,hlin,
     +     hcol,hval)

      implicit none

C     PARAMETERS
      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     SCALAR ARGUMENTS
      integer np,hnnz
      double precision lambmin,lambmax,d

C     ARRAY ARGUMENTS
      integer hlin(*),hcol(*)
      double precision s(np),x(2*np),tr(np),dtr(2*np),hval(*)

C     LOCAL SCALARS
      integer i
      double precision eta,eta2,kappa,kappa2,s2,lambi,phi,xx,xx2,tmp,AA,
     +     B1,B2,BB,C1,C2,C3,C4,C5,C6,CC,D1,D2,DD,Ax,beta,dAAde,dB1de,
     +     dB2de,dBBde,dC1de,dC2de,dC4de,dC5de,dC6de,dCCde,dD1de,dD2de,
     +     dDDde,dAxde,dbetade,dAAdk,dB1dk,dB2dk,dBBdk,dC1dk,dC2dk,
     +     dC3dk,dC5dk,dC6dk,dCCdk,dD1dk,dD2dk,dDDdk,dAxdk,dbetadk,dxdk,
     +     dAAdede,dBBdede,dC4dede,dC6dede,dCCdede,dDDdede,dAxdede,
     +     dbetadede,dBBdedk,dC6dedk,dCCdedk,dDDdedk,dAxdedk,dbetadedk,
     +     dAAdkdk,dBBdkdk,dC3dkdk,dCCdkdk,dDDdkdk,dAxdkdk,dbetadkdk

      hnnz = 0

C     Compute T[lamb_i,s(lamb_i),d,eta(lamb_i),kappa(lamb_i)] = 
C             T[lamb_i,s_i,d,eta_i,kappa_i] = tr
 
      do i = 1, np
         eta    = x(i)
         eta2   = eta**2
         kappa  = x(i+np)
         kappa2 = kappa**2
         s2     = s(i)**2
         
         lambi  = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
         phi    = (4.0d0 * pi * d * eta)/lambi
         xx     = exp((-4.0d0 * pi * kappa * d)/lambi)
         xx2    = xx**2
         
         AA     = 16.0d0 * s(i) * (eta2 + kappa2)
         
         B1     = (eta + 1.0d0)**2 + kappa2
         B2     = (eta + 1.0d0) * (eta + s2) + kappa2
         BB     = B1 * B2
         
         C1     = (eta2 + kappa2 - 1.0d0)
         C2     = (eta2 + kappa2 - s2)
         C3     = 2.0d0 * kappa2 * (s2 + 1.0d0)
         C4     = 2.0d0 * cos(phi)
         C5     = (s2 + 1.0d0) * C1
         C6     = 2.0d0 * sin(phi) * kappa
         CC     = (C1 * C2 - C3) * C4 - (2.0d0 * C2 + C5) * C6
         
         D1     = (eta - 1.0d0)**2 + kappa2
         D2     = (eta - 1.0d0) * (eta - s2) + kappa2
         DD     = D1 * D2
         
         tr(i)  = (AA * xx) / (BB - CC * xx + DD * xx2)
         
         Ax     = AA * xx
         beta   = (BB - CC * xx + DD * xx2)
         
C        First derivatives with respect to eta
            
         dAAde   = 32.0d0 * s(i) * eta
         
         dB1de   = 2.0d0 * (eta + 1.0d0)
         dB2de   = 2.0d0 * eta + s2 + 1.0d0
         dBBde   = B1 * dB2de + dB1de * B2
         
         dC1de   = 2.0d0 * eta
         dC2de   = 2.0d0 * eta
         dC4de   = - (8.0d0 * sin(phi) * pi * d )/lambi
         dC5de   = 4.0d0 * eta
         dC6de   = (8.0d0 * kappa * cos(phi) * d * pi)/lambi
         dCCde   = (C1 * C2 - C3) * dC4de + (C1*dC2de + dC1de*C2) * C4 
     +           - (2.0d0*C2 + C5) * dC6de - (2.0d0*dC2de + dC5de) * C6
         
         dD1de   = 2.0d0 * (eta - 1.0d0)
         dD2de   = 2.0d0 * eta - s2 - 1.0d0
         dDDde   = D1 * dD2de + dD1de * D2
         
         dAxde   = dAAde * xx
         dbetade = dBBde - dCCde * xx + dDDde * xx2
         
         dtr(i)  = (dAxde * beta - Ax * dbetade) / beta**2
         
C        First derivatives with respect to kappa

         dAAdk   = 32.0d0 * s(i) * kappa
         
         dB1dk   = 2.0d0 * kappa
         dB2dk   = 2.0d0 * kappa
         dBBdk   = B1 * dB2dk + dB1dk * B2
         
         dC1dk   = 2.0d0 * kappa
         dC2dk   = 2.0d0 * kappa
         dC3dk   = 4.0d0 * kappa * (s2 + 1.0d0)
         dC5dk   = 4.0d0 * kappa
         dC6dk   = 2.0d0 * sin(phi)
         dCCdk   = (C1 * dC2dk + dC1dk * C2 - dC3dk) * C4 
     +           - (2.0d0*C2 + C5) * dC6dk - (2.0d0*dC2dk + dC5dk) * C6
         
         dD1dk   = 2.0d0 * kappa
         dD2dk   = 2.0d0 * kappa
         dDDdk   = D1 * dD2dk + dD1dk * D2
         
         dxdk    = (-4.0d0 * pi * d * xx)/lambi
         
         dAxdk   = dAAdk * xx + AA * dxdk
         dbetadk = dBBdk - (dCCdk*xx + CC*dxdk) 
     +           + (dDDdk*xx2 + DD*2.0d0*xx*dxdk)
         
         dtr(i+np) = (dAxdk * beta - Ax * dbetadk) / beta**2
         
C        Second derivatives with respect to eta and eta

         dAAdede   = 32.0d0 * s(i)
         
         dBBdede   = 2.0d0 * (dB1de * dB2de + B1 + B2)
         
         dC4dede   = - (32.0d0 * cos(phi) * pi**2 * d**2 )/lambi**2
         dC6dede   = -(32.0d0*kappa*sin(phi)*d**2*pi**2)/lambi**2
         dCCdede   = (dC1de * C2 + C1 * dC2de) * dC4de 
     +             + (C1*C2 - C3) * dC4dede 
     +             + 2.0d0*C4*(dC1de * dC2de + C1 + C2) 
     +             + dC4de * (C1 * dC2de + dC1de * C2) 
     +             - 2.0d0 * dC6de * (2.0d0*dC2de+dC5de) 
     +             - (2.0d0*C2 + C5) * dC6dede - 8.0d0 * C6
         
         dDDdede   = 2.0d0 * (dD1de * dD2de + D1 + D2)
         
         dAxdede   = dAAdede * xx 
         dbetadede = dBBdede - dCCdede*xx + dDDdede*xx2
            
         tmp = (dAxdede*beta - Ax*dbetadede) * beta**2 
     +       - 2.0d0*beta*dbetade * (dAxde*beta - Ax*dbetade)

         hnnz = hnnz + 1
         hlin(hnnz) = i
         hcol(hnnz) = i
         hval(hnnz) = tmp/(beta**4)
         
C        Second derivatives with respect to eta and kappa

         dBBdedk   = dB1dk * dB2de + dB1de * dB2dk
         
         dC6dedk   = (8.0d0 * cos(phi) * d * pi)/lambi
         dCCdedk   = (dC1dk * C2 + C1 * dC2dk - dC3dk) * dC4de 
     +             + (dC1dk * dC2de + dC1de * dC2dk) * C4 
     +             - (2.0d0 * dC2dk + dC5dk) * dC6de 
     +             - (2.0d0 * C2 + C5) * dC6dedk 
     +             - (2.0d0 * dC2de + dC5de) * dC6dk
         
         dDDdedk   = dD1dk * dD2de + dD1de * dD2dk
         
         dAxdedk   = dAAde * dxdk
         dbetadedk = dBBdedk - dCCdedk*xx - dCCde*dxdk + dDDdedk*xx2 
     +             + 2.0d0*dDDde*xx*dxdk
         
         tmp = (dAxdedk*beta + dAxde*dbetadk - dAxdk*dbetade 
     +        - Ax*dbetadedk) * beta**2 - 2.0d0*beta*dbetadk 
     +        * (dAxde*beta - Ax*dbetade)
         
         hnnz = hnnz + 1
         hlin(hnnz) = i+np
         hcol(hnnz) = i
         hval(hnnz) = tmp/(beta**4)
         
C        Second derivatives with respect to kappa and kappa
         
         dAAdkdk   = 32.0d0 * s(i)
         
         dBBdkdk   = 2.0d0 * (dB1dk*dB2dk + B1 + B2)
         
         dC3dkdk   = 4.0d0 * (s2 + 1.0d0)
         dCCdkdk   = 2.0d0*C4*(dC1dk * dC2dk + C1 + C2) - dC3dkdk * C4 
     +             - 2.0d0 * dC6dk * (2.0d0*dC2dk + dC5dk) - 8.0d0 * C6
         
         dDDdkdk   = 2.0d0 * (dD1dk*dD2dk + D1 + D2)
         
         dAxdkdk   = dAAdkdk * xx + 2.0d0 * dAAdk * dxdk
         dbetadkdk = dBBdkdk - dCCdkdk*xx - 2.0d0*dCCdk*dxdk+dDDdkdk*xx2 
     +             + 2.0d0*dDDdk*xx*dxdk + 2.0d0*dxdk*(dDDdk*xx+DD*dxdk)
         
         tmp = (dAxdkdk*beta - Ax*dbetadkdk) * beta**2 
     +       - 2.0d0*beta*dbetadk * (dAxdk*beta - Ax*dbetadk)

         hnnz = hnnz + 1
         hlin(hnnz) = i+np
         hcol(hnnz) = i+np
         hval(hnnz) = tmp/(beta**4)
      end do
      
      end

C     ******************************************************************
C     ******************************************************************

      subroutine geninp(np,lambmin,lambmax,e1,e2,x)

C     SCALAR ARGUMENTS
      integer np
      double precision lambmin, lambmax, e1, e2

C     ARRAY ARGUMENTS
      double precision x(2*np)

C     LOCAL SCALARS
      integer i, j
      double precision lambi, lambk, a, b

C     Define initial values of eta and kappa

      x(np+1) = 0.1d+00
      x(2*np) = 1.0d-10
      lambk   = lambmin + 0.2d0 * (lambmax - lambmin)

      a = -0.45d0/(lambmax - lambmin)
      b = 0.1d0 - lambmin * a
      i = 2

 10   continue

      lambi = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
      if ((lambi .lt. lambk) .and. (i .le. np-1)) then
         x(np+i) = a * lambi + b
         i = i + 1
         go to 10
      end if

      a = (1.0d-02 - 1.0d-10)/(0.8d0 * lambmin - 0.8d0 * lambmax)
      b = 1.0d-10 - lambmax * a

      do j = i, np-1
         lambi = lambmin + (j-1)*(lambmax-lambmin)/(np-1)
         x(np+j) = a * lambi + b
      end do

C     Possible values of (e1,e2) = (3,2), (4,2), (5,2), (4,3), (5,3), (5,4)

      do i = 1, np
         x(i) = e1 - (e1 - e2)*(i-1) / (np-1)
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine genprob(dkick)

C     PARAMETERS
      integer    npmax
      parameter (npmax = 200)

      double precision pi
      parameter       (pi = 3.141592653589793116d0)

C     COMMON SCALARS
      integer indlaminfl,np
      double precision lambmin,lambmax,d,trued

C     COMMON ARRAYS
      double precision tobs(npmax),s(npmax),talpha(npmax),tkappa(npmax),
     +     teta(npmax)

C     LOCAL SCALARS
      integer i
      double precision lambi, lambi2, alphai, ei

C     LOCAL ARRAYS
      double precision x(2*npmax)

C     SCALAR ARGUMENTS
      double precision dkick

C     COMMON BLOCKS
      common /sprobdata/ d,trued,lambmin,lambmax,indlaminfl,np
      save   /sprobdata/
      common /vprobdata/ tobs,s,talpha,tkappa,teta
      save   /vprobdata/

      np      = 100
      lambmin =  5.5d2
      lambmax = 15.3d2
      trued   =  1.0d2

C------------------------------------------------------
C     Compute all values
C------------------------------------------------------

C     Define dkick

      dkick = trued

C     Define true values of eta, kappa and s

      do i = 1, np
         lambi    = lambmin + (i-1)*(lambmax-lambmin)/(np-1)
         lambi2   = lambi**2

         ei       = 1.24d3/lambi

         if ((ei .gt. 0.6d0) .and. (ei .le. 1.4d0)) then
            alphai = exp(6.5944d-6 * exp(9.0846d0*ei) - 16.102d0)
         elseif ((ei .gt. 1.4d0) .and. (ei .le. 1.75d0)) then
            alphai = exp(2.0d1 * ei - 41.9d0)
         elseif ((ei .gt. 1.75d0) .and. (ei .lt. 2.29d0)) then
            alphai = exp(sqrt(59.56d0 * ei - 102.1d0) - 8.391d0)
         else
            write(*,*) 'Ei: ', ei, ' undefined'
         end if

         talpha(i) = alphai
         tkappa(i) = (lambi * alphai) / (4.0d0 * pi)
         teta(i)   = sqrt(1.0d0 + 1.0d0/(9.195d-2 - 12.6d3/lambi2))
         s(i)      = sqrt(1.0d0 + 1.0d0/(7.568d-1 - 7.93d3/lambi2))
         x(i)      = teta(i)
         x(i+np)   = tkappa(i)
      end do

C     Compute true values of T
      call compt(np,lambmin,lambmax,s,trued,x,tobs)

C------------------------------------------------------
C     Write transmition values
C------------------------------------------------------

      open(11, file="prob.dat")

      do i = 1, np
         write(11,100) tobs(i)
      end do

      close(11)

 100  format(1X, F20.4)

      end
