SUBROUTINE STLIO2(ELAM,LEFT,RIGHT,MATCHP,SIZE,TOL,BOUND_ROUT,PROB_DEF,YSOL,ZSOL,LSOL,CN) ! This routine computes the eigenfunctions of 2th order sturm liouville ! problems, with the conjugate solution and the condition number. ! ! N.B. : Y is the solution of the second order problem : ! L(lambda)Y = 0 ! and Z is the solution of the problem : ! L*(lambda)Z = 0 ! ! PARAMETERS : ! ----------- ! ! *- ELAM : Complex(dble precision) ! On entry, ELAM is the eigenvalue that the user wants to use ! for the computaion. ELAM is unchanged on exit. ! *- LEFT : Real(dble precision) ! On entry, LEFT is the left interval end of integration. ! LEFT is unchanged on exit. ! *- RIGHT : Real(dble precision) ! On entry, RIGHT is the right interval end of integration. ! RIGHT is unchanged on exit. ! *- MATCHP : Real(dble precision) ! On entry, MATCHP is the matching point the user wants to ! define for the integration problem. ! MATCHP is unchanged on exit. ! Constraint : LEFT < MATCHP <= RIGHT ! If MATCHP=RIGHT, the STLIO2 routine doesn't use a matching ! point. ! *- SIZE : Integer ! On entry, SIZE is approximately the double of the size of ! solution array that the user wants to obtain. This size will ! be adjust by the routine during the computation. ! If the value supplied by the user is not big enough, a failure ! message permit to warn the user, who can try with a bigger value. ! SIZE is unchanged on exit. ! *** Recommended Value: 10000 *** ! *- TOL : Real(dble precision) ! On entry, TOL is the scalar tolerance define by the user ! for the integration. ! TOL is unchanged on exit. ! Constraint : TOL > 10.0 * machine precision ! *- BOUND_ROUT : Subroutine, supplied by the user. ! BOUND_ROUT gives the values of y and of its first derivative ! at interval ends. ! Its specification is : !----------------------------------------------------------------------- ! BOUND ROUT(IEND,X,ELAM,Y,PDY) ! ! *- IEND : integer (= 0 or 1) ! On entry, IEND indicates if the BC are needed at the left ! interval end (IEND=0), or at the right interval end (IEND=1) ! IEND is unchanged on exit. ! *- X : Real (Dble precision) ! *- ELAM : Complex (Dble precision) ! On entry, the eigenvalue of the problem, which must be the ! same as the value given to STLIO2. ! ELAM is unchanged on exit. ! *- Y : Complex (Dble precision) ! On exit, Y is the value of the solution the user wants to ! compute at the 'IENDth' end of the interval. ! *- PDY : Complex (Dble precision) ! On exit, PDY is the value of the derivative of the solution ! the user wants to compute at the 'IENDth' end of the interval. !------------------------------------------------------------------------ ! *- PROB_DEF : Subroutine, supplied by the user. ! PROB_DEF gives the values of the differential equation ! coefficients P, Q3, Q2, Q1, Q0 which defines the problem ! for one value of X between LEFT and RIGHT ! Its specification is : !------------------------------------------------------------------------ ! PROB_DEF(X,ELAM,P,Q1,Q0) ! ! *- X : Real (Dble precision) ! On entry, X is the value at which the coeffients of the ! differential equation are needed. ! *- ELAM : Complex (Dble precision) ! On entry, the eigenvalue of the problem, which must be the ! same as the value given to STLIO2. ! ELAM is unchanged on exit. ! *- P : Real (Dble precision) ! *- Q1,Q0 : Complex (Dble precision) ! On exit, P, Q3, Q2, Q1, Q0 are the differential equation ! coefficients which define the problem. !------------------------------------------------------------------------ ! *- YSOL : Real (double precision) array (0:2,1:SIZE) ! On exit, YSOL contains the computation results for Y with in its first ! row the different times at which a solution has been computed, in its ! second row the real part of the solution at this different times, and ! the last row the imaginary part of the solution at this times. ! *- ZSOL : Real (double precision) array (0:2,1:SIZE) ! On exit, ZSOL contains the computation results for Z with in its first ! row the different times at which a solution has been computed, in its ! second row the real part of the solution at this different times, and ! the last row the imaginary part of the solution at this times. ! *- LSOL : Integer ! On exit, LSOL is the useful size of the solution array. ! LSOL is approximately equal to SIZE/2. ! *- CN : Real (Dble precision) ! On exit, CN is the condition number of our problem. ( CN >= 1.0 ) IMPLICIT NONE !..... Intent arguments ......! Integer, intent(in) :: SIZE Complex(kind=kind(1.d0)), intent(in) :: ELAM Real(kind=kind(1.d0)), intent(in) :: LEFT, RIGHT, MATCHP, TOL Real(kind=kind(1.d0)), intent(out) :: CN Integer, intent(out) :: LSOL Real(kind=kind(1.d0)), Dimension(0:2,1:SIZE), intent(out) :: YSOL, ZSOL EXTERNAL BOUND_ROUT, PROB_DEF !..... Integer ......! Integer :: LA, LB, SIZEP !...... Real ........! Real(kind=kind(1.d0)) :: eps Real(kind=kind(1.d0)), Dimension(:,:), allocatable :: YSOLA, YSOLB, ZSOLA, ZSOLB Real(kind=kind(1.d0)), Dimension(:,:), allocatable :: YPSOL, YPSOLA, YPSOLB, ZPSOL, ZPSOLA, ZPSOLB !..... Complex ....... Complex(kind=kind(1.d0)) :: Y,PDY !.....Interface......! INTERFACE SUBROUTINE CALCEIGF(PI,PF,TOL,SIZE,YSOL,YPSOL,ZSOL,ZPSOL,A,LAMBDA,k,BOUND_ROUT,PROB_DEF) Integer :: NEQ, NEQG, LATOL, LRTOL, IFAIL, IEND, LRWORK, LIWORK, MAXSTP Integer, dimension(:), allocatable :: IWORK Real(kind=kind(1.d0)) :: TCRIT, TNOW, TOL, HMAX, TOUT, TINC Real(kind=kind(1.d0)), dimension(:), allocatable :: ATOL, RTOL, YNOW, YPNOW, RWORK Logical :: VECTOL, ONESTP, CRIT, ALTERG, SOPHST Character(1) :: STATEF Integer, intent(in) :: SIZE Complex(kind=kind(1.d0)) :: Y, PDY Complex(kind=kind(1.d0)), intent(in) :: LAMBDA Real(kind=kind(1.d0)), intent(in) :: PI, PF, A Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(out) :: YSOL,ZSOL,YPSOL,ZPSOL Integer, intent(out) :: k Real(kind=kind(1.d0)) :: eps, var, ceiling Integer :: i,j EXTERNAL PROB_DEF, BOUND_ROUT, F END SUBROUTINE SUBROUTINE F(T,ELAM,PROB_DEF,Y,YP) Real(kind=kind(1.d0)), intent(in) :: T Real(kind=kind(1.d0)) :: P Complex(kind=kind(1.d0)) :: Q1, Q0 Complex(kind=kind(1.d0)), intent(in) :: ELAM Real(kind=kind(1.d0)),dimension(:), intent(in) :: Y Real(kind=kind(1.d0)),dimension(:), intent(out) :: YP EXTERNAL PROB_DEF END SUBROUTINE SUBROUTINE ROUT_INTEG(F,ELAM,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL) Integer, intent(in) :: NEQ, LRWORK, LIWORK Real(kind=kind(1.d0)), dimension(LRWORK), intent(inout) :: RWORK Integer, dimension(LIWORK), intent(inout) :: IWORK Integer, intent(out) :: IFAIL Real(kind=kind(1.d0)), intent(inout) :: TNOW, TOUT Real(kind=kind(1.d0)), dimension(NEQ), intent(inout) :: YNOW Complex(kind=kind(1.d0)), intent(in) :: ELAM EXTERNAL PROB_DEF, F Integer :: NEQG, IREVCM, YRVCM, YPRVCM, KGRVCM, i Logical :: ROOT Real(kind=kind(1.d0)) :: TRVCM, GRVCM Real(kind=kind(1.d0)), dimension(NEQ) :: Y, YP END SUBROUTINE SUBROUTINE ADJUST(SIZE,YSOLA,YPSOLA,LA,YSOLB,YPSOLB,LB,YSOL,YPSOL,L) Integer, intent(in) :: SIZE Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(out) :: YSOL,YPSOL Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(in) :: YSOLA,YSOLB, YPSOLA,YPSOLB Integer, intent(in) :: LA, LB Integer,intent(out) :: L Integer :: i Real(kind=kind(1.d0)) :: norme,normep,normuse,normusep Complex(kind=kind(1.d0)) :: beta Complex(kind=kind(1.d0)),dimension(1:SIZE) :: YACOMP, YBCOMP, YPACOMP, YPBCOMP Logical :: info END SUBROUTINE SUBROUTINE CONDNUM(SIZE,YSOL,ZSOL,L,CN) Integer, intent(in) :: SIZE Real(kind=kind(1.d0)), intent(out) :: CN Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(in) :: YSOL,ZSOL Integer, intent(in) :: L Integer :: i Real(kind=kind(1.d0)) :: SY, SZ, SYZ, SYZR, SYZI Real(kind=kind(1.d0)),dimension(1:SIZE) :: TIME,YBY, ZBZ, YBZR, YBZI END SUBROUTINE SUBROUTINE INTEG(SIZE,TIME,TAB,RES,L) Integer, intent(in) :: SIZE Real(kind=kind(1.d0)), intent(out) :: RES Real(kind=kind(1.d0)),dimension(1:SIZE), intent(in) :: TIME, TAB Integer, intent(in) :: L END SUBROUTINE INTEG END INTERFACE !.... PROGRAM ........ PROGRAM ........ PROGRAM ........ PROGRAM ........ PROGRAM .....! eps = 0.001 !...... Allocate arrays .......! SIZEP=SIZE/2 Allocate(YSOLA(0:2,1:SIZEP)) Allocate(YSOLB(0:2,1:SIZEP)) Allocate(ZSOLA(0:2,1:SIZEP)) Allocate(ZSOLB(0:2,1:SIZEP)) Allocate(YPSOL(0:2,1:SIZEP)) Allocate(YPSOLA(0:2,1:SIZEP)) Allocate(YPSOLB(0:2,1:SIZEP)) Allocate(ZPSOL(0:2,1:SIZEP)) Allocate(ZPSOLA(0:2,1:SIZEP)) Allocate(ZPSOLB(0:2,1:SIZEP)) !....... Executable statements ...........! Call CALCEIGF(LEFT,MATCHP,TOL,SIZEP,YSOLA,YPSOLA,ZSOLA,ZPSOLA,LEFT,ELAM,LA,BOUND_ROUT,PROB_DEF) LSOL=LA IF (abs(MATCHP-RIGHT)<=eps) THEN YSOL = YSOLA ZSOL = ZSOLA YPSOL = YPSOLA ZPSOL = ZPSOLA ELSE Call CALCEIGF(RIGHT,MATCHP,TOL,SIZEP,YSOLB,YPSOLB,ZSOLB,ZPSOLB,LEFT,ELAM,LB,BOUND_ROUT,PROB_DEF) Call ADJUST(SIZEP,YSOLA,YPSOLA,LA,YSOLB,YPSOLB,LB,YSOL,YPSOL,LSOL) Call ADJUST(SIZEP,ZSOLA,ZPSOLA,LA,ZSOLB,ZPSOLB,LB,ZSOL,ZPSOL,LSOL) END IF Call CONDNUM(SIZE,YSOL,ZSOL,LSOL,CN) END SUBROUTINE !--------------------------------------------------------------------------------- SUBROUTINE CALCEIGF(PI,PF,TOL,SIZE,YSOL,YPSOL,ZSOL,ZPSOL,A,LAMBDA,k,BOUND_ROUT,PROB_DEF) ! This routine computes eigenfunctions ! ! N.B. : Y is the solution of the second order problem : ! L(lambda)Y = 0 ! and Z is the solution of the problem : ! L*(lambda)Z = 0 ! ! ! PARAMETERS : ! ------------ ! *- PI : real (input) : on entry, PI is the point where starts ! the integration. (unchanged on exit) ! *- PF : real (output) : on entry, PF is the point where ends ! the integration. (unchanged on exit) ! *- TOL : real (input) : on entry, TOL is the error tolerance we have ! to define for the NAG routines. ! (unchanged on exit) ! *- YSOL(0:2,1:SIZE) (output) : real array. ! This YSOL table contains in its first row the different values at ! which a solution has been computed, in its second row an approximation ! to the real part of the solution Y at each time of the first row, and ! in its third row an approximation to the imaginary part of the solution ! Y at the same time. ! *- YPSOL(0:2,1:SIZE) (output) : real array. ! This YPSOL table contains in its first row the different values at ! which a solution has been computed, in its second row an approximation ! to the derivate of the real part of the solution Y at each time of the ! first row, and in its third row an approximation to the derivate of ! the imaginary part of the solution Y at the same time. ! *- ZSOL(0:2,1:SIZE) (output) : real array. ! This ZSOL table contains in its first row the different values at ! which a solution has been computed, in its second row an approximation ! to the real part of the solution Z at each time of the first row, and ! in its third row an approximation to the imaginary part of the solution ! Z at the same time. ! *- ZPSOL(0:2,1:SIZE) (output) : real array. ! This ZPSOL table contains in its first row the different values at ! which a solution has been computed, in its second row an approximation ! to the derivate of the real part of the solution Z at each time of the ! first row, and in its third row an approximation to the derivate of ! the imaginary part of the solution Z at the same time. ! *- A : real (input) : on entry, A is a reference point,the left interval end ! which permit to define the problem (boundary conditions, ! the sense of the integration). (unchanged on exit) ! *- LAMBDA : real (input) : on entry, Lambda is the parameter needed to define ! the operator L. The different possible values of ! LAMBDA are the eigenvalues given by the routine. ! *- k : integer (output) : k give us a very important information, the quantity ! of datas stored in the different arrays. So, the ! parameter SIZE has to be greater than k. Otherwise, ! the routine fails with an error message to warn the ! user. ! *- SIZE : integer (input) : SIZE is the SOL arrays size chosen by the user before ! to call this routine. ! WARNING : it is not the quantity of datas in output. ! *- BOUND_ROUT : external routine (input) : BOUND_ROUT is the routine supplied ! by the user which defines the boudary conditions at left ! and right interval ends. ! *- PROB_DEF : external routine (input) : PROB_DEF is the routine supplied ! by the user which defines the differential equation ! coefficients. IMPLICIT NONE !.......... Variables for the NAG routines .........! Integer :: NEQ, NEQG, LATOL, LRTOL, IFAIL, IEND, LRWORK, LIWORK, MAXSTP Integer, dimension(:), allocatable :: IWORK Real(kind=kind(1.d0)) :: TCRIT, TNOW, HMAX, TOUT, TINC Real(kind=kind(1.d0)), dimension(:), allocatable :: ATOL, RTOL, YNOW, YPNOW, RWORK Logical :: VECTOL, ONESTP, CRIT, ALTERG, SOPHST Character(1) :: STATEF Real(kind=kind(1.d0)), intent(in) :: TOL !........ Scalars & Arrays ..........! Integer, intent(in) :: SIZE Complex(kind=kind(1.d0)) :: Y, PDY Complex(kind=kind(1.d0)), intent(in) :: LAMBDA Real(kind=kind(1.d0)), intent(in) :: PI, PF, A Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(out) :: YSOL,ZSOL,YPSOL,ZPSOL Integer, intent(out) :: k Real(kind=kind(1.d0)) :: eps, var, ceiling Integer :: i,j EXTERNAL PROB_DEF, BOUND_ROUT, F !.... Parameters....! NEQ = 8 eps = TOL ceiling = 10.0 STATEF = 'S' LATOL = 1 LRTOL = 1 VECTOL = .false. ONESTP = .false. CRIT = .true. TCRIT = PF HMAX = 0.0 MAXSTP = 0 NEQG = 0 ALTERG = .false. SOPHST = .false. LRWORK = 23*(1+NEQ)+4 LIWORK = 21 + NEQG TINC = 2.D0*(PF-PI)/SIZE !.... Arrays .....! Allocate(YNOW(NEQ)) Allocate(YPNOW(NEQ)) Allocate(ATOL(LATOL)) Allocate(RTOL(LRTOL)) Allocate(RWORK(LRWORK)) Allocate(IWORK(LIWORK)) ATOL(1) = TOL RTOL(1) = TOL !............ Executable statements .................! IF ((PI-A) <= eps) THEN IEND = 0 ELSE IEND = 1 END IF Call BOUND_ROUT(IEND,PI,LAMBDA,Y,PDY) YNOW(1) = Dble(Y) YNOW(2) = Dimag(Y) YNOW(3) = Dble(PDY) YNOW(4) = Dimag(PDY) YNOW(5) = Dble(Y) YNOW(6) = -Dimag(Y) YNOW(7) = Dble(PDY) YNOW(8) = -Dimag(PDY) IFAIL = 0 Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT,HMAX,MAXSTP,NEQG,& & ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL) TNOW = PI k=1 YSOL(0,k) = TNOW YSOL(1,k) = YNOW(1) YSOL(2,k) = YNOW(2) YPSOL(0,k) = TNOW YPSOL(1,k) = YNOW(3) YPSOL(2,k) = YNOW(4) ZSOL(0,k) = TNOW ZSOL(1,k) = YNOW(5) ZSOL(2,k) = YNOW(6) ZPSOL(0,k) = TNOW ZPSOL(1,k) = YNOW(7) ZPSOL(2,k) = YNOW(8) DO WHILE (abs(TNOW-PF)>eps) IF (k >= SIZE) THEN Write(*,*) 'The table dimension is not big enough, try with a bigger value.' STOP ELSE IFAIL = -1 TOUT = TNOW + TINC IF ( PI < PF .and. TOUT > TCRIT) THEN TOUT = PF ELSE IF ( PI > PF .and. TOUT < TCRIT) THEN TOUT = PF END IF Call ROUT_INTEG(F,LAMBDA,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL) IF (IFAIL==0) THEN IF ((PI-A) <= eps) THEN YSOL(0,k+1) = TNOW YSOL(1,k+1) = YNOW(1) YSOL(2,k+1) = YNOW(2) YPSOL(0,k+1) = TNOW YPSOL(1,k+1) = YNOW(3) YPSOL(2,k+1) = YNOW(4) ZSOL(0,k+1) = TNOW ZSOL(1,k+1) = YNOW(5) ZSOL(2,k+1) = YNOW(6) ZPSOL(0,k+1) = TNOW ZPSOL(1,k+1) = YNOW(7) ZPSOL(2,k+1) = YNOW(8) ELSE DO i=1,k j=k-i+1 YSOL(0,j+1) = YSOL(0,j) YSOL(1,j+1) = YSOL(1,j) YSOL(2,j+1) = YSOL(2,j) YPSOL(0,j+1) = YPSOL(0,j) YPSOL(1,j+1) = YPSOL(1,j) YPSOL(2,j+1) = YPSOL(2,j) ZSOL(0,j+1) = ZSOL(0,j) ZSOL(1,j+1) = ZSOL(1,j) ZSOL(2,j+1) = ZSOL(2,j) ZPSOL(0,j+1) = ZPSOL(0,j) ZPSOL(1,j+1) = ZPSOL(1,j) ZPSOL(2,j+1) = ZPSOL(2,j) END DO YSOL(0,1) = TNOW YSOL(1,1) = YNOW(1) YSOL(2,1) = YNOW(2) YPSOL(0,1) = TNOW YPSOL(1,1) = YNOW(3) YPSOL(2,1) = YNOW(4) ZSOL(0,1) = TNOW ZSOL(1,1) = YNOW(5) ZSOL(2,1) = YNOW(6) ZPSOL(0,1) = TNOW ZPSOL(1,1) = YNOW(7) ZPSOL(2,1) = YNOW(8) END IF var = max(abs(YNOW(1)),abs(YNOW(2)),abs(YNOW(3)),abs(YNOW(4)),abs(YNOW(5)),abs(YNOW(6)),abs(YNOW(7)),abs(YNOW(8))) IF (var >= ceiling) THEN DO i=1,k+1 YSOL(1,i) = YSOL(1,i)/var YSOL(2,i) = YSOL(2,i)/var YPSOL(1,i) = YPSOL(1,i)/var YPSOL(2,i) = YPSOL(2,i)/var ZSOL(1,i) = ZSOL(1,i)/var ZSOL(2,i) = ZSOL(2,i)/var ZPSOL(1,i) = ZPSOL(1,i)/var ZPSOL(2,i) = ZPSOL(2,i)/var END DO YNOW(1) = YNOW(1)/var YNOW(2) = YNOW(2)/var YNOW(3) = YNOW(3)/var YNOW(4) = YNOW(4)/var YNOW(5) = YNOW(5)/var YNOW(6) = YNOW(6)/var YNOW(7) = YNOW(7)/var YNOW(8) = YNOW(8)/var IF (abs(TNOW-PF)>=eps) THEN STATEF = 'R' Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT, & & HMAX,MAXSTP,NEQG,ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL) END IF ELSE STATEF = 'C' Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT, & & HMAX,MAXSTP,NEQG,ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL) END IF ELSE Write(*,*) 'IFAIL =', IFAIL END IF k=k+1 END IF END DO END SUBROUTINE CALCEIGF !-------------------------------------------------------------------------------- SUBROUTINE ROUT_INTEG(F,ELAM,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL) IMPLICIT NONE !...... Intent arguments ........! Integer, intent(in) :: NEQ, LRWORK, LIWORK Real(kind=kind(1.d0)), dimension(LRWORK), intent(inout) :: RWORK Integer, dimension(LIWORK), intent(inout) :: IWORK Integer, intent(out) :: IFAIL Real(kind=kind(1.d0)), intent(inout) :: TNOW, TOUT Real(kind=kind(1.d0)), dimension(NEQ), intent(inout) :: YNOW Complex(kind=kind(1.d0)), intent(in) :: ELAM EXTERNAL PROB_DEF, F !....... Local arguments ..........! Integer :: NEQG, IREVCM, YRVCM, YPRVCM, KGRVCM, i Logical :: ROOT Real(kind=kind(1.d0)) :: TRVCM, GRVCM Real(kind=kind(1.d0)), dimension(NEQ) :: Y, YP !........ Executable statements ........! NEQG = 0 ROOT = .false. IREVCM = 0 IFAIL = -1 Call D02QGF(NEQ,TNOW,YNOW,TOUT,NEQG,ROOT,IREVCM,TRVCM,YRVCM,YPRVCM,GRVCM,KGRVCM,RWORK,LRWORK,IWORK,LIWORK,IFAIL) DO WHILE (IREVCM > 0 .and. IREVCM < 9) IF( IREVCM < 8 ) THEN IF (YRVCM == 0) THEN Y = YNOW ELSE DO i=1,NEQ Y(i) = RWORK(YRVCM+i-1) END DO END IF Call F(TRVCM,ELAM,PROB_DEF,Y,YP) DO i=1,NEQ RWORK(YPRVCM+i-1) = YP(i) END DO END IF IFAIL = -1 Call D02QGF(NEQ,TNOW,YNOW,TOUT,NEQG,ROOT,IREVCM,TRVCM, & & YRVCM,YPRVCM,GRVCM,KGRVCM,RWORK,LRWORK,IWORK,LIWORK,IFAIL) END DO END SUBROUTINE ROUT_INTEG !-------------------------------------------------------------------------------- SUBROUTINE F(T,ELAM,PROB_DEF,Y,YP) ! This routine is the routine needed by the NAG routine ROUT_INTEG to evaluate y prime IMPLICIT NONE !.... Scalar arguments .....! Real(kind=kind(1.d0)), intent(in) :: T Real(kind=kind(1.d0)) :: P Complex(kind=kind(1.d0)), intent(in) :: ELAM Complex(kind=kind(1.d0)) :: Q1, Q0 !.... Array arguments ......! Real(kind=kind(1.d0)),dimension(8), intent(in) :: Y Real(kind=kind(1.d0)),dimension(8), intent(out) :: YP !....... Other arguments .........! EXTERNAL PROB_DEF !........ Executable statements .........! Call PROB_DEF(T,ELAM,P,Q1,Q0) YP(1) = Y(3)/P YP(2) = Y(4)/P YP(3) = ((Dble(Q1))*(Y(3))-(Dimag(Q1))*(Y(4)))/P+(Dble(Q0))*Y(1)-(Dimag(Q0))*Y(2) YP(4) = ((Dble(Q1))*(Y(4))+(Dimag(Q1))*(Y(3)))/P+(Dble(Q0))*Y(2)+(Dimag(Q0))*Y(1) YP(5) = Y(7)/P YP(6) = Y(8)/P YP(7) = ((-Dble(Q1))*(Y(7))-(Dimag(Q1))*(Y(8)))/P+(Dble(Q0))*Y(5)+(Dimag(Q0))*Y(6) YP(8) = ((-Dble(Q1))*(Y(8))+(Dimag(Q1))*(Y(7)))/P+(Dble(Q0))*Y(6)-(Dimag(Q0))*Y(5) END SUBROUTINE F !---------------------------------------------------------------------------------------- SUBROUTINE ADJUST(SIZE,YSOLA,YPSOLA,LA,YSOLB,YPSOLB,LB,YSOL,YPSOL,L) IMPLICIT NONE !....... Local arrays ........! Integer, intent(in) :: SIZE Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(out) :: YSOL,YPSOL Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(in) :: YSOLA,YSOLB, YPSOLA,YPSOLB !...... Local integers .......! Integer, intent(in) :: LA, LB Integer,intent(out) :: L Integer :: i !........ Real arguments ........! Real(kind=kind(1.d0)) :: norme,normep,normuse,normusep !........ Complex arguments & arrays ........! Complex(kind=kind(1.d0)) :: beta Complex(kind=kind(1.d0)),dimension(1:SIZE) :: YACOMP, YBCOMP, YPACOMP, YPBCOMP !....... Logical argument...........! Logical :: info !........ Executable statements ........! norme = ((YSOLA(1,LA)**2+YSOLB(1,1)**2)+(YSOLA(2,LA)**2+YSOLB(2,1)**2))**(0.5) normep = ((YPSOLA(1,LA)**2+YPSOLB(1,1)**2)+(YPSOLA(2,LA)**2+YPSOLB(2,1)**2))**(0.5) DO i=1,LA YACOMP(i) = DCMPLX(YSOLA(1,i),YSOLA(2,i)) END DO DO i=1,LB YBCOMP(i) = DCMPLX(YSOLB(1,i),YSOLB(2,i)) END DO DO i=1,LA YPACOMP(i) = DCMPLX(YPSOLA(1,i),YPSOLA(2,i)) END DO DO i=1,LB YPBCOMP(i) = DCMPLX(YPSOLB(1,i),YPSOLB(2,i)) END DO IF ( norme >= normep ) THEN normuse = (YSOLA(1,LA)**2+YSOLA(2,LA)**2)**(0.5) normusep = (YSOLB(1,1)**2+YSOLB(2,1)**2)**(0.5) IF( normuse >= normusep ) THEN beta = YBCOMP(1)/YACOMP(LA) info = .true. ELSE beta = YACOMP(LA)/YBCOMP(1) info = .false. END IF ELSE normuse = (YPSOLA(1,LA)**2+YPSOLA(2,LA)**2)**(0.5) normusep = (YPSOLB(1,1)**2+YPSOLB(2,1)**2)**(0.5) IF( normuse >= normusep ) THEN beta = YPBCOMP(1)/YPACOMP(LA) info = .true. ELSE beta = YPACOMP(LA)/YPBCOMP(1) info = .false. END IF END IF IF (info) THEN DO i=1,LA YACOMP(i)=YACOMP(i)*beta YPACOMP(i)=YPACOMP(i)*beta END DO ELSE DO i=1,LB YBCOMP(i)=YBCOMP(i)*beta YPBCOMP(i)=YPBCOMP(i)*beta END DO END IF DO i=1,LA YSOL(0,i) = YSOLA(0,i) YSOL(1,i) = Dble(YACOMP(i)) YSOL(2,i) = Dimag(YACOMP(i)) YPSOL(0,i) = YPSOLA(0,i) YPSOL(1,i) = Dble(YPACOMP(i)) YPSOL(2,i) = Dimag(YPACOMP(i)) END DO DO i=2,LB YSOL(0,LA+i-1) = YSOLB(0,i) YSOL(1,LA+i-1) = Dble(YBCOMP(i)) YSOL(2,LA+i-1) = Dimag(YBCOMP(i)) YPSOL(0,LA+i-1) = YPSOLB(0,i) YPSOL(1,LA+i-1) = Dble(YPBCOMP(i)) YPSOL(2,LA+i-1) = Dimag(YPBCOMP(i)) END DO L=LA+LB-1 END SUBROUTINE ADJUST !----------------------------------------------------------------------------------------- SUBROUTINE CONDNUM(SIZE,YSOL,ZSOL,L,CN) ! This routine calculate the condition number using the INTEG routine. ! The function Y and Z are define by the entry tables YSOL and ZSOL, ! and the result is the real number CN (on exit). IMPLICIT NONE !........... Intent arguments ..............! Integer, intent(in) :: SIZE Real(kind=kind(1.d0)), intent(out) :: CN Real(kind=kind(1.d0)),dimension(0:2,1:SIZE), intent(in) :: YSOL,ZSOL Integer, intent(in) :: L !........... Local arguments ...........! Integer :: i Real(kind=kind(1.d0)) :: SY, SZ, SYZR, SYZI, SYZ Real(kind=kind(1.d0)),dimension(1:SIZE) :: TIME,YBY, ZBZ,YBZR,YBZI !.......... Executable statements .............! DO i=1,L TIME(i) = YSOL(0,i) YBY(i) = YSOL(1,i)**2+YSOL(2,i)**2 ZBZ(i) = ZSOL(1,i)**2+ZSOL(2,i)**2 YBZR(i) = YSOL(1,i)*ZSOL(1,i)+YSOL(2,i)*ZSOL(2,i) YBZI(i) = YSOL(1,i)*ZSOL(2,i)-YSOL(2,i)*ZSOL(1,i) END DO Call INTEG(SIZE,TIME,YBY,SY,L) Call INTEG(SIZE,TIME,ZBZ,SZ,L) Call INTEG(SIZE,TIME,YBZR,SYZR,L) Call INTEG(SIZE,TIME,YBZI,SYZI,L) SYZ = (SYZR**2+SYZI**2)**(0.5) CN = (SY**(0.5)*SZ**(0.5))/SYZ END SUBROUTINE CONDNUM !---------------------------------------------------------------------------------------- SUBROUTINE INTEG(SIZE,TIME,TAB,RES,L) ! This routine integrate the entry function defined by TIME and TAB, using the ! trapezoidal rule, and RES is the the result of the integration, on exit. IMPLICIT NONE !........... Intent arguments ..............! Integer, intent(in) :: SIZE Real(kind=kind(1.d0)), intent(out) :: RES Real(kind=kind(1.d0)),dimension(1:SIZE), intent(in) :: TIME, TAB Integer, intent(in) :: L !........... Local arguments ...........! Integer :: i !.......... Executable statements .............! RES = 0.0 DO i=1,L-1 RES = RES+(0.5*(TAB(i+1)+TAB(i))*(TIME(i+1)-TIME(i))) END DO END SUBROUTINE INTEG !------------------------------------------------------------------------------------