stirling.f

Go to the documentation of this file.

!                       *******************
                        SUBROUTINE stirling
!                       *******************
!
     &(ni,xi,yi,no,xostep,yo)
!
!***********************************************************************
! ARTEMIS   V7P0                                   June 2014
!***********************************************************************
!
!brief    INTERPOLATES A FUNCTION TO DESIRED RESOLUTION.
!
!history  N.DURAND (HRW)
!+        March 2001
!+
!+   Ported to V7P0
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| NI XI YI       |   | INPUT NUMBER OF DISCRETE VALUES,
!|                |   | INPUT X AND Y VALUES
!| NO XOSTEP YO   |   | OUTPUT NUMBER OF DISCRETE VALUES,
!|                |   | CORRESPONDING OUTPUT RESOLUTION IN X,
!|                |   | AND OUTPUT Y VALUES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          :: NI,NO
      DOUBLE PRECISION :: XI(ni),YI(ni)
      DOUBLE PRECISION :: XOSTEP,XO,YO(no)
!
      INTEGER          :: I,IV,IVP
      DOUBLE PRECISION :: A1,A2,A3,A4,A5,A6
      DOUBLE PRECISION :: RMIN,R,R2,U
!
      INTRINSIC float, dabs
!
!-----------------------------------------------------------------------
      xostep = (xi(ni) - xi(1)) / float(no-1)
      yo(1)  = yi(1)
      yo(no) = yi(ni)
!
!-----------------------------------------------------------------------
!     PARAMETERS FOR THE PARABOLA FITTING
!     OUTSIDE DO AND IF LOOPS FOR OPTIMISATION PURPOSES
!     ESTIMATED ONLY ONCE
      a3 = ( yi(3)-yi(1) )/( xi(3)-xi(1) )/( xi(3)-xi(2) ) -
     &     ( yi(2)-yi(1) )/( xi(2)-xi(1) )/( xi(3)-xi(2) )
      a2 = ( yi(2)-yi(1) )/( xi(2)-xi(1) ) - a3*( xi(2)+xi(1) )
      a1 = yi(1) - ( a2 + a3*xi(1) )*xi(1)
!
      a4 = ( yi(ni-2)-yi(ni) )/( xi(ni-2)-xi(ni) )/( xi(ni-2)-xi(ni-1) )
     &   - ( yi(ni-1)-yi(ni) )/( xi(ni-1)-xi(ni) )/( xi(ni-2)-xi(ni-1) )
      a5 = ( yi(ni-1)-yi(ni) )/( xi(ni-1)-xi(ni) )
     &   - a4*( xi(ni-1)+xi(ni) )
      a6 = yi(ni) - ( a5 + a4*xi(ni) )*xi(ni)
!
!-----------------------------------------------------------------------
!     LOOP ON THE INTERPOLATION POINTS IN BETWEEN
      ivp = 2
      DO i = 2,(no-1)
        rmin = 1.d0
        xo = xi(1)+float(i-1)*xostep
!
!-----------------------------------------------------------------------
!       DETERMINES IVP SUCH THAT XO IS CLOSEST TO XI(IVP)
!       R<0 FOR XO BEFORE XI(IVP); R>0 FOR XO AFTER XI(IVP)
        DO iv = ivp,(ni-1)
          r = ( xo - xi(iv) ) / ( xi(iv+1) - xi(iv) )
          IF( dabs(r).LE.dabs(rmin) ) THEN
            rmin = r
            ivp = iv
          ENDIF
        ENDDO ! IV = IVP,(NI-1)
        r = ( xo - xi(ni) ) / ( xi(ni) - xi(ni-1) )
        IF( dabs(r).LE.dabs(rmin) ) THEN
          rmin = r
          ivp = ni
        ENDIF
        r = rmin
!
!-----------------------------------------------------------------------
!       XO CLOSEST TO XI(2) :
!       YO COMPUTED FROM PARABOLA USING 1ST 3 POINTS
        IF( ivp.EQ.2 ) THEN
          yo(i) = a1 + a2*xo + a3*xo**2
!
!-----------------------------------------------------------------------
!       XO CLOSEST TO XI(3) OR XI(NI-2) :
!       2ND ORDER STIRLING WITH CONSTANT STEP
        ELSEIF( ivp.EQ.3 .OR. ivp.EQ.ni-2 ) THEN
          r2 = r**2
          yo(i) = ( r2-r )*yi(ivp-1)/2.d0 +
     &        ( 1.d0-r2 )*yi(ivp)  +  ( r2+r )*yi(ivp+1)/2.d0
!
!-----------------------------------------------------------------------
!       XO CLOSEST TO XI(NI-1) OR XI(NI) :
!       YO COMPUTED FROM PARABOLA USING LAST 3 POINTS
        ELSEIF( ivp.EQ.ni-1 .OR. ivp.EQ.ni ) THEN
          yo(i) = a6 + a5*xo + a4*xo**2
!
!-----------------------------------------------------------------------
!       SUBSEQUENT POINTS
!       4TH ORDER STIRLING WITH CONSTANT STEP
        ELSE
          u = r*( (r**2)-1.d0 )/12.d0
          yo(i) =     u*( (r/2.d0)-1.d0 ) * yi(ivp-2) +
     &        (1.d0-r)*(2.d0*u-(r/2.d0))  * yi(ivp-1) +
     &          ( r*(3.d0*u-r) + 1.d0 )   * yi(ivp)   +
     &        (1.d0+r)*(-2.d0*u+(r/2.d0)) * yi(ivp+1) +
     &                u*( (r/2.d0)+1.d0 ) * yi(ivp+2)
        ENDIF
!
!-----------------------------------------------------------------------
!     END OF THE LOOP ON THE INTERPOLATION POINTS
      ENDDO ! I = 2,(NO-1)
!
      RETURN
      END