NAG Library Manual, Mark 30.1
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   G13EJF Example Program Text
!   Mark 30.1 Release. NAG Copyright 2024.

    Module g13ejfe_mod

!     G13EJF Example Program Module:
!     User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: f, h, read_problem_data
!     .. Parameters ..
      Integer, Parameter, Public       :: mx = 3, my = 2, nin = 5, nout = 6
!     .. Derived Type Definitions ..
      Type, Public                     :: g13ej_problem_data
        Real (Kind=nag_wp)             :: delta, a, r, d
        Real (Kind=nag_wp)             :: phi_rt, phi_lt
      End Type g13ej_problem_data
    Contains
      Subroutine f(n,xt,fxt,dat)

!       .. Scalar Arguments ..
        Type (g13ej_problem_data), Intent (In) :: dat
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: fxt(:,:)
        Real (Kind=nag_wp), Intent (In) :: xt(:,:)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: t1, t3
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, sin
!       .. Executable Statements ..

        t1 = 0.5_nag_wp*dat%r*(dat%phi_rt+dat%phi_lt)
        t3 = (dat%r/dat%d)*(dat%phi_rt-dat%phi_lt)

        Do i = 1, n
          fxt(1,i) = xt(1,i) + cos(xt(3,i))*t1
          fxt(2,i) = xt(2,i) + sin(xt(3,i))*t1
          fxt(3,i) = xt(3,i) + t3
        End Do

        Return
      End Subroutine f
      Subroutine h(n,yt,hyt,dat)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Type (g13ej_problem_data), Intent (In) :: dat
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: hyt(:,:)
        Real (Kind=nag_wp), Intent (In) :: yt(:,:)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tmp
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, sin
!       .. Executable Statements ..

        Do i = 1, n
          hyt(1,i) = dat%delta - yt(1,i)*cos(dat%a) - yt(2,i)*sin(dat%a)
          hyt(2,i) = yt(3,i) - dat%a

!         Make sure that the theta is in the same range as the observed
!         data, which in this case is [0, 2*pi)
          If (hyt(2,i)<0.0_nag_wp) Then
            hyt(2,i) = hyt(2,i) + 2*x01aaf(tmp)
          End If
        End Do

        Return
      End Subroutine h
      Subroutine read_problem_data(t,dat,read_ok)
!       Read in any data specific to the F and H subroutines

!       .. Scalar Arguments ..
        Type (g13ej_problem_data), Intent (Inout) :: dat
        Integer, Intent (In)           :: t
        Logical, Intent (Out)          :: read_ok
!       .. Local Scalars ..
        Integer                        :: tt
!       .. Executable Statements ..

        If (t==0) Then
!         Read in the data that is constant across all time points
          Read (nin,*) dat%r, dat%d, dat%delta, dat%a
          read_ok = .True.
        Else
!         Read in data for time point t
          Read (nin,*) tt, dat%phi_rt, dat%phi_lt
          If (tt/=t) Then
!           Sanity check
            Write (nout,99999) 'Expected to read in data for time point ', t
            Write (nout,99999) 'Data that was read in was for time point ', tt
99999       Format (A,E22.15)
            read_ok = .False.
          Else
            read_ok = .True.
          End If
        End If
      End Subroutine read_problem_data
    End Module g13ejfe_mod

    Program g13ejfe

!     .. Use Statements ..
      Use g13ejfe_mod, Only: f, g13ej_problem_data, h, mx, my, nin, nout,      &
                             read_problem_data
      Use nag_library, Only: g13ejf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (g13ej_problem_data)        :: dat
      Integer                          :: i, ifail, irevcm, ldfxt, ldlx, ldly, &
                                          ldst, ldxt, licomm, lrcomm, lropt,   &
                                          n, ntime, t
      Logical                          :: read_ok
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: fxt(:,:), lx(:,:), ly(:,:),          &
                                          rcomm(:), ropt(:), st(:,:), x(:),    &
                                          xt(:,:), y(:)
      Integer, Allocatable             :: icomm(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, max, repeat
!     .. Executable Statements ..
      Write (nout,*) 'G13EJF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)

!     Using default optional arguments
      lropt = 0
      Allocate (ropt(lropt))

!     Allocate arrays
      n = 2*mx + 1
      If (lropt>=1) Then
        If (abs(ropt(1)-2.0_nag_wp)<=0.0_nag_wp) Then
          n = n + 2*mx
        End If
      End If
      ldlx = mx
      ldly = my
      ldst = mx
      ldxt = mx
      ldfxt = max(mx,my)
      licomm = 30
      lrcomm = 30 + my + mx*my + 2*max(mx,my)
      Allocate (lx(ldlx,mx),ly(ldly,my),x(mx),st(ldst,mx),xt(ldxt,max(my,      &
        n)),fxt(ldfxt,n+max(mx,my)),icomm(licomm),rcomm(lrcomm),y(my))

!     Read in the Cholesky factorization of the covariance matrix for the
!     process noise
      Do i = 1, mx
        Read (nin,*) lx(i,1:i)
      End Do

!     Read in the Cholesky factorization of the covariance matrix for the
!     observation noise
      Do i = 1, my
        Read (nin,*) ly(i,1:i)
      End Do

!     Read in the initial state vector
      Read (nin,*) x(1:mx)

!     Read in the Cholesky factorization of the initial state covariance
!     matrix
      Do i = 1, mx
        Read (nin,*) st(i,1:i)
      End Do

!     Read in the number of time points to run the system for
      Read (nin,*) ntime

!     Read in any problem specific data that is constant
      Call read_problem_data(0,dat,read_ok)
      If (.Not. read_ok) Then
        Go To 100
      End If

!     Title for first set of output
      Write (nout,*) ' Time   ', repeat(' ',(11*mx-16)/2), 'Estimate of State'
      Write (nout,*) repeat('-',7+11*mx)

!     Loop over each time point
      irevcm = 0
      Do t = 1, ntime

!       Read in any problem specific data that is time dependent
        Call read_problem_data(t,dat,read_ok)
        If (.Not. read_ok) Then
          Go To 100
        End If

!       Read in the observed data for time t
        Read (nin,*) y(1:my)

!       Call Unscented Kalman Filter routine
ukf_lp: Do
          ifail = 0
          Call g13ejf(irevcm,mx,my,y,lx,ldlx,ly,ldly,x,st,ldst,n,xt,ldxt,fxt,  &
            ldfxt,ropt,lropt,icomm,licomm,rcomm,lrcomm,ifail)
          Select Case (irevcm)
          Case (1)
!           Evaluate F(X)
            Call f(n,xt,fxt,dat)

          Case (2)
!           Evaluate H(X)
            Call h(n,xt,fxt,dat)

          Case Default
!           IREVCM = 3, finished
            Exit ukf_lp
          End Select
        End Do ukf_lp

!       Display the current state estimates
        Write (nout,99999) t, x(1:mx)
      End Do

      Write (nout,*)
      Write (nout,*) 'Estimate of Cholesky Factorization of the State'
      Write (nout,*) 'Covariance Matrix at the Last Time Point'
      Do i = 1, mx
        Write (nout,99998) st(i,1:i)
      End Do

100   Continue

99999 Format (1X,I3,4X,10(1X,F10.3))
99998 Format (10(1X,E10.3))
    End Program g13ejfe