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