c                                        
c	do_MRK_calc  
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Controlling routine that calls the marker calculation
c				routines
c
c	Routines called:
c		InterpolateMrkData 
c		AdjustMrkData
c		SmoothMrkData
c		DigFilterMrkData
c		RegularizeMrkData
c     MrkPosSegLCS
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param			GP.analysis_type
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	None
c
      subroutine do_MRK_calc( CMND,MF,SF,JF,CF,OF )
      implicit none

      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_names.inc'

      record /CmndLine/     CMND
      record /marker_flag/  MF                
      record /segment_flag/ SF
      record /joint_flag/   JF
      record /cycle_flag/   CF
      record /option_flag/  OF

      if ( CMND.CmndText(2) .eq. 'INT' ) then 
        call InterpolateMrkData( CMND,MF,OF )

      else if ( CMND.CmndText(2) .eq. 'ADJ' ) then
        if ( (GP.analysis_type .eq. FULL_ANL) .or.
     &       (GP.analysis_type .eq. HH_ANL)       ) then
          call AdjustMrkData( MF, OF )
        else if ( GP.analysis_type .eq. CKC_ANL ) then
          call AdjustCKCMrkData( MF, OF )
        else
          write(6,110) ANLNames(GP.analysis_type)
          if ( OF.ToLog ) write(3,110) ANLNames(GP.analysis_type)
110       format(' Error- Marker adjustment is not needed for ',a20,' analysis')
        end if                              

      else if ( CMND.CmndText(2) .eq. 'SMO' ) then
        call SmoothMrkData( CMND,MF,OF )

      else if ( CMND.CmndText(2) .eq. 'FIL' ) then
        call DigFilterMrkData( CMND,MF,OF )

      else if ( CMND.CmndText(2) .eq. 'REG' ) then
        call RegularizeMrkData( CMND,MF,SF,JF,CF,OF )

      else if ( CMND.CmndText(2) .eq. 'SLP' ) then
        call MrkPosSegLCS( CMND,MF,SF,OF )

      else                                  
        write(6,100) CMND.CmndText(2)
        if ( OF.ToLog ) write(3,100) CMND.CmndText(2)
100     format(' ',a3,' is an invalid MRK command.'/
     &         ' Valid commands are: INTerpolate ADJust SMOoth'/
     &         '                     REGularize  FILter'/
     &         '                     SegmentLocalPosition')
      end if      
      return
      end                  
c                                        
c	InterpolateMrkData  
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Find marker positions which are obviously bad due to 
c				drop out more than probably from the motion measurement
c				system.  Having found these interpolate the missing
c				positions back using a polynomial of selectable order.
c				The routine will also extrapolate marker positions, 
c				using a different polynomial if desired, with due 
c				warning that the resulting marker positions are more 
c				than probably fiction rather than fact.  When this 
c				occurs one must be very careful evaluating any subsequent
c				quantities calculated at or near these extrapolated 
c				frames which are dependent upon the marker position
c				data.  A good way of looking at the error is to 
c				calculate the segment intermarker distances and see if
c				vary alot at these frames (SEG MRK).  NOTE that this is a
c				sufficient but NOT necessary condition for trouble.  One
c				last comment: Each file format/motion system has a 
c				different way of indicating a bad frame for a marker
c				and hence this routine works for certain with only 
c				TR3 and C3D files since these are the only ones that I 
c				have had access to so far. The function bad_marker is
c				the key to finding bad marker positions so look there
c				if you are going to support a new file format...  PS,
c				please send me the code too :-)
c
c	Routines called:
c		bad_marker 
c		find_interp_frm
c		lagrange_interp
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param			GP.analysis_type
c	Param			GP.int_poly_order
c	Param			rMtnFrmRate
c	Param			MP.mrk_name
c	Param			MP.nfrm
c	Param			MP.num_mrk
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	Marker			MrkData
c
      subroutine InterpolateMrkData( CMND,MF,OF )
      implicit none

      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_mrk.inc'
      include 'a_names.inc'
                              
      record /CmndLine/     CMND
      record /marker_flag/  MF
      record /option_flag/  OF
                                                             
      integer*4 poly_order,mrk,num_bad_frm,num_good_frm,frm,bad_frm(MAX_MTN),
     &          good_frm(MAX_MTN),i,interp_frm(10),axis,j,k,extra_poly_order,
     &          temp_poly_order
      real*4    time_data(10),pos_data(10,3),interp_time,interp_val
      logical*4 bad_marker,error,skip_mrk,UseExtraPoly
                                              
        
      MF.interpolated= .false.
      poly_order= GP.int_poly_order ! maximum of 9
      extra_poly_order= poly_order
      UseExtraPoly= .false.
c                                                                           
c     Check for temporary override of polynomial orders- both interpolation and
c     extrapolation
c
      if ( CMND.NumOpt(2) .gt. 0 ) then
        do i= 1, CMND.NumOpt(2) 
          if ( CMND.OptText(i,2) .eq. 'POL' ) then
            if ( (.not. CMND.OptNumOK(i,2)) .or.
     &           ( CMND.OptIntVal(i,2) .lt. 1) .or.
     &           ( CMND.OptIntVal(i,2) .gt. 9) ) goto 1010
            poly_order= CMND.OptIntVal(i,2)
            extra_poly_order= poly_order
          else if ( CMND.OptText(i,2) .eq. 'EXT' ) then
            if ( (.not. CMND.OptNumOK(i,2)) .or.
     &           ( CMND.OptIntVal(i,2) .lt. 1) .or.                
     &           ( CMND.OptIntVal(i,2) .gt. 9) ) goto 1020
            temp_poly_order= CMND.OptIntVal(i,2)
            UseExtraPoly= .true.
          else                                 
            goto 1030
          end if                          
        end do ! i                             
      end if ! options are present
c
c     First check for bad marker frames
c
      if ( UseExtraPoly ) extra_poly_order= temp_poly_order ! this is cumbersome
      if ( MF.read ) then                              ! but assures that 
        do mrk= 1, MP.num_mrk         ! the extra poly will be assigned
          skip_mrk= .false.                ! no matter the order of the options
          do i= 1, num_init_empty_mrk      
            if ( mrk .eq. init_empty_mrk(i) ) skip_mrk= .true.
          end do
          if ( .not. skip_mrk ) then                  
            if ( OF.UpDate ) write(6,101) mrk,MP.num_mrk
101         format('+Currently checking marker # ',i2,'/',i2,'.')
            num_bad_frm= 0
            num_good_frm= 0
*$ifdef DEBUG
            write(3,444) mrk
 444        format(' Marker: ',i3)
*$endif
            do frm= 1, MP.nfrm
              if ( bad_marker(MrkData(x,frm,mrk)) ) then
                num_bad_frm= num_bad_frm + 1
                bad_frm(num_bad_frm)= frm 
              else
                num_good_frm= num_good_frm + 1
                good_frm(num_good_frm)= frm
              end if    
            end do ! frm
            if ( num_bad_frm .gt. 0 ) then
              if ( num_bad_frm .le. 4 ) then
                if ( OF.UpDate ) then
                  write(6,100) mrk,MP.num_mrk,MP.mrk_name(mrk),
     &                         (bad_frm(i),i=1,num_bad_frm)
                  if ( OF.ToLog ) write(3,100) mrk,MP.num_mrk,
     &                                               MP.mrk_name(mrk),
     &                                               (bad_frm(i),i=1,num_bad_frm)
100              format('+Currently checking marker # ',i2,'/',i2,
     &                 '  Bad data in ',a4,' frame(s): ',5i4/)
                end if
              else
                if ( num_bad_frm .le. 11 ) then
                  write(6,105) MP.mrk_name(mrk),(bad_frm(i),i=1,num_bad_frm)
                  if ( OF.ToLog ) write(3,105) MP.mrk_name(mrk),
     &                                                 (bad_frm(i),i=1,num_bad_frm)
105               format(' Bad data in ',a4,' frame(s): ',11i4/)
                else
                  write(6,102) MP.mrk_name(mrk),(bad_frm(i),i=1,11)
                  if ( OF.ToLog ) write(3,102) MP.mrk_name(mrk),
     &                                                 (bad_frm(i),i=1,11)
102               format(' Bad data in ',a4,' frame(s): ',11i4)
                  j= 12
                  do while ( j .le. num_bad_frm )
                    k= j + 11
                    if ( k .gt. num_bad_frm ) k= num_bad_frm
                    write(6,103) (bad_frm(i),i=j,k)
                    if ( OF.ToLog ) write(3,103) (bad_frm(i),i=j,k)
103                 format(28x,12i4)
                    j= j + 12
                  end do
                  write(6,104)
                  if ( OF.ToLog ) write(3,104)
104               format(/)                  
                end if ! num_bad_frm gt 11
              end if ! num_bad_frm gt 4
            end if ! are there any bad frames ?     
c
c           if there are bad frames then interpolate/extrapolate them
c
            if ( (num_bad_frm .gt. 0) .and. 
     &           (num_good_frm .lt. poly_order+1) ) then
              write(6,200) mrk
              if ( OF.ToLog ) write(3,200) mrk
200           format(' ERROR- Marker #',i2,' can be interpolated because'/
     &               '        there are not enough good frames of data'/
     &               '        The marker data will be left as is...'/)
            else
              do frm= 1, num_bad_frm
c
c      Check whether extrapolation will occur and if so switch to extra_poly
c
                if ( (bad_frm(frm) .lt. good_frm(1)) .or. 
     &               (bad_frm(frm) .gt. good_frm(num_good_frm)) ) then
                  temp_poly_order= extra_poly_order
                  write(6,300) bad_frm(frm)
                  if ( OF.ToLog ) write(3,300) bad_frm(frm)
300               format(' WARNING- Extrapolating frame # ',i4)
                else
                  temp_poly_order= poly_order
                end if
                call find_interp_frm( good_frm,num_good_frm,bad_frm(frm),
     &                                temp_poly_order+1,interp_frm )
                do i= 1, temp_poly_order+1
                  time_data(i)= float(interp_frm(i))/rMtnFrmRate
                  do axis= x, z
                    pos_data(i,axis)= MrkData(axis,interp_frm(i),mrk)
                  end do ! axis
                end do ! i            
                interp_time= float(bad_frm(frm))/rMtnFrmRate
                do axis= x, z
                  call lagrange_interp( time_data,pos_data(1,axis),
     &                                  temp_poly_order+1,
     &                                  interp_time,interp_val,error )
                  if ( error ) goto 1000
                  MrkData(axis,bad_frm(frm),mrk)= interp_val     
                end do ! axis
              end do ! frm
            end if ! too many frames missing
          end if ! .not. skip_mrk
        end do ! mrk   
        MF.interpolated= .true.
        if ( GP.analysis_type .eq. UPR_ANL ) MF.adjusted= .true.
      else
        write(6,110)                   
        if ( OF.ToLog ) write(3,110)
110     format(' Marker data has not been read in yet.  Interpolation can not'/
     &         ' take place until the marker data is present.  Use the REAd'/
     &         ' command to read in marker data from a file.')
      end if
      return           
                      
1000  write(6,1005) mrk,bad_frm(frm)
      if ( OF.ToLog ) write(3,1005) mrk,bad_frm(frm)
1005  format(' Error- An error occured while interpolating marker # ',i2,
     &       ' at frame # ',i3,'.' )
      return

1010  write(6,1015) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,1015) CMND.OptText(i,2)
1015  format(' Error- ',a3,' is not a valid interpolation polynomial order'/
     &       '        Valid orders are from 1 to 9')
      return

1020  write(6,1025) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,1025) CMND.OptText(i,2)
1025  format(' Error- ',a3,' is not a valid extrapolation polynomial order'/
     &       '        Valid orders are from 1 to 9')
      return

1030  write(6,1035) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,1035) CMND.OptText(i,2)
1035  format(' Error- ',a3,' is not a valid option for MaRKer INTerpolate'/
     &       '        Valid options are : POLynomialorder (interp/extrap)'/
     &       '                            EXTrapolation polyorder')
      return
      end         
c                                        
c	bad_marker  
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Identify marker position as one containing drop out 
c				in the position data.  This routine only supports
c				TR3 (VICON) and C3D (AMASS) files at present.
c           
c	Routines called:
c
c	VecMag
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	None
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	None
c
      function bad_marker( marker )
      implicit none

      real*4    marker(3),VecMag,MrkMag
      logical*4 bad_marker

      bad_marker= .false.
      MrkMag= VecMag(marker)
*$ifdef DEBUG
      write(3,100) MrkMag
 100  format(f10.6)
*$endif
      if ( MrkMag .lt. 0.0001 ) then ! drop out in tr3 files
        bad_marker= .true. 
      else if ( (abs(marker(1)+1.0) .lt. 1.0e-5) .and.
     &          (abs(marker(2)+1.0) .lt. 1.0e-5) .and.
     &          (abs(marker(3)+1.0) .lt. 1.0e-5)       ) then
        bad_marker= .true.                   ! drop out in c3d files
      else if ( (abs(marker(1)) .lt. 1.0e-5) .and.
     &          (abs(marker(2)) .lt. 1.0e-5) .and.
     &          (abs(marker(3))-0.483 .lt. 1.0e-5) ) then
        bad_marker= .true.                   ! drop out in c3d files
      end if
      return
      end
c                                        
c	find_interp_frm  
c
c	Copyright by Dwight Meglan 1990
c
c
c	Purpose:	Find frames of marker data that can be used to 
c				interpolate or extrapolate a marker's position
c           
c	Routines called:
c
c	None
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	None
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	None
c
      subroutine find_interp_frm( good_pnt_num, num_good_pnts, bad_pnt_num,
     &                            num_int_pnts, interp_pnts )
      implicit none
                                      
      integer*4 num_good_pnts,good_pnt_num(num_good_pnts),bad_pnt_num,
     &          num_int_pnts,interp_pnts(num_int_pnts)
 
      integer*4 i,num_lt_good_pnts,num_gt_good_pnts,nearest_lower_pnt,
     &          num_lt_int_pnts,num_gt_int_pnts,lower_good_pnt_num,
     &          upper_good_pnt_num,pnt,num_pnts
c                                                      
c     Find the number of data points greater than and less than the x
c     interpolation value.                  
c
      i= 1
      do while ( (good_pnt_num(i).lt.bad_pnt_num) .and. (i.lt.num_good_pnts) )
        i= i + 1                      
      end do                
      num_lt_good_pnts= i - 1
      nearest_lower_pnt= i - 1
      num_gt_good_pnts= num_good_pnts - num_lt_good_pnts
c                  
c     Decide which data points can be used for the interpolation.  The
c     following set of if statements assumes that there are at least as
c     many data points as there are interpolation points required.  This
c     algorithm nicely takes care of the bad data points lying at the beginning
c     or the end of the data set by adjusting the interpolation point numbers
c     such that they are all either after (in the case of backward
c     extrapolation) or before (in the case of forward extrapolation) the
c     point number being extrapolated.
c
      num_lt_int_pnts= num_int_pnts/2
      num_gt_int_pnts= num_int_pnts - num_lt_int_pnts
      if ( num_lt_good_pnts .lt. num_lt_int_pnts ) then
        num_lt_int_pnts= num_lt_good_pnts                ! decrease lt_pnts
        num_gt_int_pnts= num_int_pnts - num_lt_good_pnts ! increase gt_pnts
      else if ( num_gt_good_pnts .lt. num_gt_int_pnts ) then
        num_gt_int_pnts= num_gt_good_pnts                ! decrease gt_pnts
        num_lt_int_pnts= num_int_pnts - num_gt_good_pnts ! increase lt_pnts
      end if     
      lower_good_pnt_num= nearest_lower_pnt - num_lt_int_pnts + 1
      upper_good_pnt_num= lower_good_pnt_num + num_int_pnts - 1
c                          
c     Assign points to be used in the interpolation
c                 
      num_pnts= 1
      do pnt= lower_good_pnt_num, upper_good_pnt_num
        interp_pnts(num_pnts)= good_pnt_num(pnt)
        num_pnts= num_pnts + 1
      end do     
      return
      end
c                                        
c	AdjustMrkData  
c
c	Copyright by Dwight Meglan 1990
c
c
c	Purpose:	This is a bit of a kludgey (sp?) routine to adjust the
c				position of some markers in the OSU full body marker set
c				and create a couple of new positions from several of
c				marker's positions.  Specificially, the knee and ankle
c				markers are moved medially by a distance based on the
c				ratio of knee/ankle width of calf length.  The hip positions
c				are created using the right and left ASIS markers and
c				the right and left knee markers.  This method is documented
c				by Tylkowski (~1983) as well as Brand (1990).  I have
c				evaluated the phantom translation at the hip due the
c				estimation of the hip location using the JNT TRN
c				command (look there for more details, its quite useful)
c				and have found it to be a maximum of 2 cm which, though
c				not perfect, is adequate I believe.  Many of the original
c				concepts come from the old Boston Children's Hospital
c				Gait Lab software with a number of modifications by moi.
c				This process is absolutely necessary for the subsequent
c				calculations to be accurate.
c           
c	Routines called:
c
c	MinLineXPnts
c	Tricen
c	VecMag
c	VecNorm
c	VecSub
c	VecX
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param			MRK_NUM
c	Param			MP.nfrm
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	Marker			MrkData
c
      subroutine AdjustMrkData( MF, OF )
      implicit none
      
      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_mrk.inc'
      include 'a_names.inc'

      record /marker_flag/  MF
      record /option_flag/  OF

      integer*4 side,frm,axis,VecNorm,MinLineXPnts,error
      real*4    sign,vecA(3),vecB(3),vecC(3),vecD(3),vecE(3),normD(3),normE(3),
     &          ankle_pnt(3),foot_pnt(3),knee_width,asis_line(3),z_axis(3),
     &          x_axis(3),thigh_length,shldr_width,shldr_line(3),asis_width,
     &          VecMag,calc_error
c     
c     Transfer the knee markers medially such that they will be situated at
c     the 'middle' of the knee.  The markers are
c     translated along a line formed from the ankle, knee and tibial stick
c     markers which is approximately medial (rt) or lateral (lf).  The distance
c     along this line that the marker is translated is proportional to the
c     distance between the knee and ankle markers.  The ratios of joint
c     width to tibial length are measured from several subjects and a linear
c     relation is assumed between these lengths.  These ratios and the linearity
c     assumption is used without proof.  For the tibial markers, a more
c     complicated method is used because the plane formed by the foot markers
c     (toe,lateral mallelous,and heel) is not parallel to the plane formed
c     by the ankle, tibial stick, and knee markers when a person is standing
c     prone.  In fact, the foot plane is such that a static pronation appears
c     to be present in the ankle unless the position of the ankle,toe, and
c     heel markers are corrected.  The method used here is explained in
c     the derivations which accompany ANALYZE.
c
      calc_error= 0
      MF.adjusted= .false.
      if ( .not. MF.interpolated ) goto 1000 ! error
      do frm= 1, MP.nfrm         
        if ( OF.UpDate .and. (mod(frm,10) .eq. 0) )  
     &    write(6,100) frm,MP.nfrm
100     format('+Processing frame # ',i3,'/',i3,'.')      
        sign= 1.0    
        do side= rt, lf
          if ( side .eq. lf ) sign= -1.0
          call VecSub( MrkData(x,frm,MRK_NUM(side).knee),
     &                  MrkData(x,frm,MRK_NUM(side).ankle),
     &                  vecA                                      )
          call VecSub( MrkData(x,frm,MRK_NUM(side).calf),
     &                  MrkData(x,frm,MRK_NUM(side).ankle),
     &                  vecB                                      )
          call VecX( vecA, vecB, vecC )
          call VecX( vecA, vecC, vecD )     
          calc_error= VecNorm( vecD, normD )
          if ( calc_error .ne. 0 ) goto 1010
          call VecSub( MrkData(x,frm,MRK_NUM(side).toe),
     &                  MrkData(x,frm,MRK_NUM(side).heel),
     &                  vecE                                     )
          calc_error= VecNorm( vecE, normE )
          if ( calc_error .ne. 0 ) goto 1010
          calc_error= MinLineXPnts( normD,MrkData(x,frm,MRK_NUM(side).ankle),
     &                              normE,MrkData(x,frm,MRK_NUM(side).heel),
     &                              ankle_pnt,foot_pnt  )
          if ( calc_error .ne. 0 ) goto 1010
          knee_width=  (5.25/42.0)*VecMag( vecA ) ! approx. 1/2 knee width
          do axis= x, z
            MrkData(axis,frm,MRK_NUM(side).knee)=
     &         MrkData(axis,frm,MRK_NUM(side).knee) + knee_width*normD(axis)
            MrkData(axis,frm,MRK_NUM(side).ankle)= ankle_pnt(axis)
          end do ! axis
c
c     Calculate the approximate location of the hip using the method developed
c     by Tylkowski and Simon
c               
          if ( side .eq. rt ) then ! do only once per frm
            call VecSub( MrkData(x,frm,MRK_NUM(lf).asis),
     &                    MrkData(x,frm,MRK_NUM(rt).asis),
     &                    vecA                                   )
            calc_error= VecNorm( vecA, asis_line )
            if ( calc_error .ne. 0 ) goto 1010
            asis_width= VecMag( vecA )
            call VecSub( MrkData(x,frm,MRK_NUM(rt).sacrum),
     &                    MrkData(x,frm,MRK_NUM(rt).asis),
     &                    vecA                                     )
            call VecX( asis_line, vecA, vecB )         
            calc_error= VecNorm( vecB, z_axis )
            if ( calc_error .ne. 0 ) goto 1010
            call VecX( asis_line, z_axis, vecB )
            calc_error= VecNorm( vecB, x_axis )
            if ( calc_error .ne. 0 ) goto 1010
          end if ! side = rt
          call VecSub( MrkData(x,frm,MRK_NUM(side).asis),
     &                  MrkData(x,frm,MRK_NUM(side).knee),
     &                  vecA                                     )
          thigh_length= VecMag( vecA )
          do axis= x, z
            MrkData(axis,frm,MRK_NUM(side).hip)=
     &                 MrkData(axis,frm,MRK_NUM(side).asis)
     &               + 0.11*sign*asis_width*asis_line(axis)
     &               - 0.21*asis_width*x_axis(axis)
     &               - 0.12*thigh_length*z_axis(axis)
          end do ! axis
        end do ! side
c
c     Set location of a neck marker at the midpoint of a line between the
c     right and left shoulder markers.
c
        call VecSub( MrkData(x,frm,MRK_NUM(lf).shldr),
     &                MrkData(x,frm,MRK_NUM(rt).shldr),
     &                vecA                                    )
        shldr_width= VecMag( vecA )
        calc_error= VecNorm( vecA, shldr_line )
        if ( calc_error .ne. 0 ) goto 1010
        do axis= x, z
          MrkData(axis,frm,MRK_NUM(rt).neck)= 
     &      MrkData(axis,frm,MRK_NUM(rt).shldr)+0.5*shldr_width*shldr_line(axis)
        end do ! axis
c
c       Calculate mid sacrum position- not used by ANZ but needed for
c        compatability with saving tru files
c
        call Tricen( MrkData(x,frm,MRK_NUM(rt).asis),
     &               MrkData(x,frm,MRK_NUM(lf).asis),
     &               MrkData(x,frm,MRK_NUM(rt).sacrum),       
     &               MrkData(x,frm,13),error             )
c        if ( error .ne. 0 ) then
c          if ( OF.UpDate ) write(6,110) frm,MP.nfrm
c110       format('+Processing frame # ',i3,'/',i3,'   Warning- Mid sacrum'/)
c        end if
      end do ! frm                  
      MF.adjusted= .true.
      return

1000  write(6,1005)
      if ( OF.ToLog ) write(3,1005)
1005  format(' Error- Marker data has not been checked and interpolated yet.'/
     &       '        Before adjustment can take place, this must be done- MRK INT.')
      return

1010  write(6,1015)
      if ( OF.ToLog ) write(3,1015)
1015  format(' Error- An error occured while calculating the ankle marker position.'/
     &       '        The marker positions were not successfully adjusted and as a'/
     &       '        result further analysis of the data is not possible.')
      return
      end      
c                                        
c	Lagrange_Interp  
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:
c		This routine performs generalized lagrangian interpolation based upon
c		the pairs of x,y data in xdata and ydata.  The interpolating polynomial
c		is of the order numpnt-1.  Note that xint, the x point at which the
c		interpolated y is calculated does not have be between the minimum
c		and maximum values of x.  Thus, in this case the y value is actually
c		an extrapolated value.  In this mode, caution must used, since extra-
c		polating too far will lead to very erroneous y's.  Lastly, note that
c		the x data points do not have to be equispaced .
c
c           
c	Routines called:
c
c	None
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	None
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	None
c
      Subroutine Lagrange_Interp( xdata, ydata, numpnt, xint, yint, error )
      implicit none
      
      integer*4 numpnt,i,j
      real*4    xdata(1),ydata(1),xint,yint,term
      logical*4 error

      error= .false.
      yint= 0.0  
      do i= 1, numpnt                                                    
        term= ydata(i)
        do j= 1, numpnt
          if ( i .ne. j ) then
            if ( abs(xdata(i)-xdata(j)) .lt. 0.00001 ) goto 1000
            term= term*( xint-xdata(j) )/( xdata(i)-xdata(j) )
          end if
        end do
        yint= yint + term
      end do
      return

1000  error= .true.
      return
      end
c                                        
c	SmoothMrkData  
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Use either the dynamic programming filter of Busby &
c				Trujillo or the Generalized Cross Validation filter
c				of Herman Woltring to smooth the marker trajectories.
c				Optionally the marker trajectory derivatives can be
c				calculated.
c           
c	Routines called:
c
c	filter
c	RemoveCMNDOpt
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param			rMtnFrmRate
c	Param			MP.nfrm
c	Param			MP.num_mrk
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	Marker			MrkAcc
c	Marker			MrkData
c	Marker			MrkVel
c
      Subroutine SmoothMrkData( CMND,MF,OF )

      implicit none
                           
      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_mrk.inc'
      include 'a_names.inc'

      record /CmndLine/     CMND        
      record /marker_flag/  MF
      record /option_flag/  OF
                            
      real*4    interval,temp(MAX_MTN),vel(MAX_MTN),acc(MAX_MTN)
      integer*4 mrk,frm,axis,status,i,j
      logical*4 SmoothOnly,DoVelAcc
      character axis_names(3)*1    
                                
      data axis_names / 'x','y','z' /

      if ( (.not. MF.adjusted) .and. 
     &     ( (GP.Analysis_Type .eq. HH_ANL) .or.
     &       (GP.Analysis_Type .eq. FULL_ANL)   ) ) goto 1000 ! error
c
c     Check for velocity/acceleration flag
c
      DoVelAcc= .false.
      SmoothOnly= .true.                
      if ( CMND.NumOpt(2) .gt. 0 ) then
        do i= CMND.NumOpt(2), 1, -1
          if ( CMND.OptText(i,2) .eq. 'VAA' ) then
            DoVelAcc= .true.
            SmoothOnly= .false.
            call RemoveCMNDOpt( CMND,2,i )
          end if ! vaa option
        end do ! i
      end if ! num opt > 0
      interval= 1.0/rMtnFrmRate
      do mrk= 1, MP.num_mrk
        if ( OF.UpDate ) write(6,100) mrk,MP.num_mrk
100     format('+Currently filtering marker # ',i2,'/',i2,'.' )
        do axis= x, z
          do frm= 1, MP.nfrm
            temp(frm)= MrkData(axis,frm,mrk)
          end do ! frm
          call filter( CMND,SmoothOnly,MP.nfrm,temp,interval,OF,
     &                 vel,acc,status )
          if ( status .lt. 0 ) goto 1010
          if ( status .eq. 4 ) then
            write(6,120) mrk,MP.num_mrk,axis_names(axis)
120         format('+Currently filtering marker # ',i2,'/',i2,
     &             '   Filter error- 'a1,' not smoothed'/)
          else if ( status .eq. 3 ) then
            write(6,110) mrk,MP.num_mrk,axis_names(axis)
110         format('+Currently filtering marker # ',i2,'/',i2,
     &             '   Warning- ',a1,' data very small'/)
          else if ( status .eq. 2 ) then
            write(6,130) mrk,MP.num_mrk,axis_names(axis)
130         format('+Currently filtering marker # ',i2,'/',i2,
     &               '   Warning- ',a1,' nearly constant'/)
          else if ( status .ge. 10 ) then
            write(6,140) mrk,MP.num_mrk,axis_names(axis),status-10
140         format('+Currently filtering marker # ',i2,'/',i2,
     &               '   Warning- ',a1,'  GCVSPL error # ',i2)
          else              
            do frm= 1, MP.nfrm
              MrkData(axis,frm,mrk)= temp(frm)
            end do ! frm                       
            if ( DoVelAcc ) then
              do frm= 1, MP.nfrm
                MrkVel(axis,frm,mrk)= vel(frm)
                MrkAcc(axis,frm,mrk)= acc(frm)
              end do
            end if ! vel/accel are to calculated
          end if            
        end do ! axis                          
      end do ! mrk
      if ( MF.NumProcOps .le. 4 ) then
        MF.NumProcOps= MF.NumProcOps + 1
        MF.ProcOp(MF.NumProcOps)= CMND.Text(1)(1:40)
      end if
      MF.smoothed= .true.
      if ( DoVelAcc ) MF.VelAcc= .true.
      return                               

1000  write(6,1005)   
      if ( OF.ToLog ) write(3,1005)
1005  format(' Error- Marker data has not been adjusted yet.'/
     &       '        Before smoothing can take place, this must be done- MRK ADJ.')
      return

1010  write(6,1015)
      if ( OF.ToLog ) write(3,1015)
1015  format(' Error - Smoothing operation unsuccessful')
      return
      
1020  write(6,1025)
      if ( OF.ToLog ) write(3,1025)
1025  format(' Error - Invalid option for marker smoothing'/
     &       '         Valid options are: /DPB /SPL')
      return                 
      end                                       
c                                        
c	DigFilterMrkData
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Provide digitial filtering of marker data if desired.
c				This is an optional facilty provided primarily for 
c				experimenting with digital filtering compared to GCV
c				 or DPF.  I don't recommend its use with clinical
c				data.  This is a controlling routine which builds the 
c				frequency domain filter function and then drives the
c				actual fft and ifft routines used for the filter.
c
c	Routines called:
c	BuildFilter
c	DigFilterDataSet
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param		   rMtnFrmRate
c	Param			MP.nfrm
c	Param			GP.num_mrk
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	Marker			MrkData
c
      Subroutine DigFilterMrkData( CMND,MF,OF )
                                                              
      implicit none
                           
      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_mrk.inc'
      include 'a_names.inc'

      record /CmndLine/     CMND        
      record /marker_flag/  MF
      record /option_flag/  OF
c
c     Zero Pad data so that the number of data points is a power of 2, apply
c     a window function to the data since it is assumed to be nonperiodic,
c     take the fourier transform of the result then multiply it with the filter
c     function, and take the inverse transform.
c
                     
      integer*4 i,j,n,NumFilPnts,NumBeginZero,NumEndZero,mrk,axis,NumBlank
      real*4    Filter(2048),smplfreq,data(MAX_MTN)
      character AxisName(3)*1
      logical*4 Error,DoWindow,BlankEnds

      data AxisName / 'x','y','z' /               
                                      
      if ( .not. MF.read ) goto 1000
      if ( rMtnFrmRate .lt. 1.0e-5 ) goto 1010 ! error
      if ( MP.nfrm .lt. 2 ) goto 1020 ! error
c
c     Figure out how many zeros must be added to data.  Zeros are added
c     symmetrically to the beginning and end of the array.  This should
c     take care of most of the end effects resulting from fourier transform
c     process.             
c                               
      n= MP.nfrm                                
      i= 1
      do while( (n - 2**i) .gt. 0 )
        i= i + 1                        
      end do                
      i= i + 1
      if ( i .gt. 11 ) goto 1030 ! error- too many samples, must be <=512
      NumFilPnts= 2**i                                                
      NumBeginZero= (NumFilPnts-n)/2
      NumEndZero=   (NumFilPnts-n) - NumBeginZero                          
c                              
c     Build filter function
c                                                       
      smplfreq= rMtnFrmRate
      call BuildFilter( CMND,2,NumFilPnts,smplfreq,
     &                  Filter,DoWindow,BlankEnds,NumBlank,Error,OF )
      if ( Error ) return       
c                   
c     Do Filtering of each marker set            
c
      do mrk= 1, MP.num_mrk
        do axis= x, z
          if ( OF.UpDate ) write(6,100) mrk,MP.num_mrk,AxisName(axis)
100       format('+ Filtering marker # ',i2,'/',i2,' ',a1)
          do i= 1, MP.nfrm
            data(i)= MrkData(axis,i,mrk)
          end do
          call DigFilterDataSet( data,n,
     &                           NumFilPnts,NumBeginZero,NumEndZero,
     &                           Filter,DoWindow,BlankEnds,NumBlank )
          do i= 1, MP.nfrm
            MrkData(axis,i,mrk)= data(i)
          end do
        end do ! axis 
      end do ! mrk               
      MF.smoothed= .true.       
      if ( MF.NumProcOps .le. 4 ) then
        MF.NumProcOps= MF.NumProcOps + 1
        MF.ProcOp(MF.NumProcOps)= CMND.Text(1)(1:40)
      end if          
      return
                 
1000  write(6,1005)
      if ( OF.ToLog ) write(3,1005)
1005  format(' Error- Marker data has not been read in yet')
      return

1010  write(6,1015) rMtnFrmRate
      if ( OF.ToLog ) write(3,1015) rMtnFrmRate
1015  format(' Error- Marker sample rate is ',f10.6,' Hz which is too small.'/
     &       '        Something has gone wrong in reading the marker file')
      return                
      
1020  write(6,1025) MP.nfrm
      if ( OF.ToLog ) write(3,1025) MP.nfrm
1025  format(' Error- There are only ',i4,' samples of data in each marker'/
     &       '        There must be at least 2 samples to perform filtering')
      return
                            
1030  write(6,1035) MP.nfrm
      if ( OF.ToLog ) write(3,1035) MP.nfrm
1035  format(' Error- There are ',i6,' samples per marker set, but at most'/
     &       '        512 samples/channel can be filtered')
      return
      end                   
c                                        
c	RegularizeMrkData
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Redistribute the markers over a new time interval and
c				at a selectable time increment.  This routine is primarily
c				intended to allow later comparision across gait cycles
c				by allowing a single gait cycle to be broken down into
c				uniform intervals which are consistent from one cycle
c				to the next.  For example, the cycle may be set so that
c				ther are 50 points over it effectively breaking it into
c				2 % of cycle intervals.  The down side of doing this is
c				that the time intervals between the points from one gait
c				cycle to next will not be the same.  Oh well, you can't
c				have both- that's physics for you.  Use this routine
c				cautiously.  SInce it merely fits an interpolating 
c				polynomial to the data, it can produce noisy results if
c				the marker data is already noisy.  Hence, I have found
c				that filter the marker data (MRK SMO), then regularizing
c 				it (MRK REG) and then filtering it again causes a minimum
c				maiming of the origianl intent of the marker trajectories.
c				Another caution is that once this routine has been run
c				all segment and joint quantities must be reevaluated
c				since the starting point data of the calculation sequence
c				has been changed.  By the way, the name regularize was
c				used to indicate that more than simple interpolation is
c				occuring, like to fill in missing data points, and that it has
c				serious implications as to affecting the marker data.
c				Lastly, note that a provision for adding 'pad' points
c				before and after the desired interpolation area is 
c				available.  This is very useful for later smooth/differentiation
c				to minimimize end effect influence of the numerical
c				techniques upon the marker data.
c
c	Routines called:
c
c	CheckGaitPattern
c	find_interp_frm
c	lagrange_interp
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c	Param			GP.int_poly_order
c	Param			GP.num_mrk
c	Param			MP.cycle_first_frm
c	Param			MP.cycle_last_frm
c	Param			MP.num_lhs
c	Param			MP.num_lto
c	Param			MP.num_rhs
c	Param			MP.num_rto
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
c	Param			rMtnFrmRate
c	Param			MP.lhs_frm
c	Param			MP.lto_frm
c	Param			MP.nfrm
c	Param			MP.rhs_frm
c	Param			MP.rto_frm
c	Marker			MrkData
c
      subroutine RegularizeMrkData( CMND,MF,SF,JF,
     &                              CF,OF )
      implicit none

      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_names.inc'
      include 'a_mrk.inc'
      include 'a_gls2.inc'
      include 'a_frc.inc'

      record /CmndLine/     CMND
      record /marker_flag/  MF
      record /segment_flag/ SF
      record /joint_flag/   JF
      record /cycle_flag/   CF
      record /option_flag/  OF

      integer*4 StartFrm,EndFrm,NumPnts,i,j,frm,mrk,GoodFrm(MAX_MTN),
     &          NumGoodFrm,NumBadFrm,NewPnt,OldFrm,PolyOrder,axis,NumPadPnts,
     &          FirstFrm,InterpFrm(10)
      real*4    NewMrkData(3,MAX_MTN,MAX_MRK),TimeData(MAX_MTN),NewFrmRate,
     &          YData(10),XData(10),FrmRateRatio,NewTime,TimeInc
      logical*4 UseCycleInfo,error,CheckGaitPattern

      if ( .not. MF.adjusted ) goto 1000 ! error
      PolyOrder= GP.int_poly_order
      UseCycleInfo= .false.
      StartFrm= 1
      EndFrm=   MP.nfrm
      NumPnts=  100
      NumPadPnts= 0
      if ( CMND.NumOpt(2) .gt. 0 ) then
        do i= 1, CMND.NumOpt(2)
          if ( CMND.OptText(i,2) .eq. 'CYC' ) then
            if ( .not. CF.event ) goto 1010 ! error
            StartFrm= MP.cycle_first_frm(1)
            EndFrm=   MP.cycle_last_frm(1)
            UseCycleInfo= .true.
          else if ( CMND.OptText(i,2) .eq. 'STA' ) then
            if ( .not. CMND.OptNumOK(i,2) ) goto 1020
            StartFrm= CMND.OptIntVal(i,2) - MP.first_frm + 1
            if ( (StartFrm .gt. MP.nfrm) .or.
     &           (StartFrm .lt. 1)                     ) goto 1020
            if ( UseCycleInfo ) goto 1060
          else if ( CMND.OptText(i,2) .eq. 'END' ) then
            if ( .not. CMND.OptNumOK(i,2) ) goto 1030
            EndFrm= CMND.OptIntVal(i,2) - MP.first_frm + 1
            if ( (EndFrm .le. StartFrm) .or.
     &           (EndFrm .gt. MP.nfrm) ) goto 1030
            if ( UseCycleInfo ) goto 1070 
          else if ( CMND.OptText(i,2) .eq. 'PNT' ) then
            if ( (.not. CMND.OptNumOK(i,2)) .or.
     &           ( CMND.OptIntVal(i,2) .lt. 1) .or.
     &           ( CMND.OptIntVal(i,2) .gt. MAX_MTN) ) goto 1040
            NumPnts= CMND.OptIntVal(i,2)
          else if ( CMND.OptText(i,2) .eq. 'POL' ) then
            if ( (.not. CMND.OptNumOK(i,2)) .or.
     &           ( CMND.OptIntVal(i,2) .lt. 1) .or.
     &           ( CMND.OptIntVal(i,2) .gt. 9) ) goto 2010
            PolyOrder= CMND.OptIntVal(i,2)
          else if ( CMND.OptText(i,2) .eq. 'PAD' ) then
            if ( (.not. CMND.OptNumOK(i,2)) .or.
     &           ( CMND.OptIntVal(i,2) .lt. 0) .or.
     &           ( CMND.OptIntVal(i,2) .gt. 10) ) goto 2020
            NumPadPnts= CMND.OptIntVal(i,2)
          else                                 
            goto 1050                     
          end if
        end do ! i
      end if ! options are present

c                                                                  
c     Perform the regularization (actually just a fancy name for interpolation)
c
      if ( rMtnFrmRate .lt. 1.0e-5 ) goto 1080    
      do i= 1, MP.nfrm
        TimeData(i)= float(i-1)/rMtnFrmRate
        GoodFrm(i)= i
      end do
      NumGoodFrm= MP.nfrm
      NumBadFrm= 0          
      NewFrmRate= rMtnFrmRate*float(NumPnts)/float(EndFrm-StartFrm+1)
      NewTime= float(StartFrm)/rMtnFrmRate
      TimeInc= 1.0/NewFrmRate
c
c     If pad points are being used to allow differentiation and smoothing later
c     on then adjust the number of points and the start time
c
      if ( NumPadPnts .gt. 0 )  NewTime= NewTime - float(NumPadPnts)*TimeInc
      do NewPnt= 1, NumPnts + NumPadPnts*2
        OldFrm= NewTime*rMtnFrmRate
        if ( OldFrm .lt. 1 ) OldFrm= 1  ! just in case NumPadPnts pushes it -
        if ( OldFrm .gt. MP.nfrm ) OldFrm= MP.nfrm
        call find_interp_frm( GoodFrm,NumGoodFrm,OldFrm,PolyOrder+1,InterpFrm )
        do j= 1, PolyOrder+1
          YData(j)= TimeData(InterpFrm(j))
        end do
c  Debug 8/1/90
c        write(10,110) NewTime,NewPnt,OldFrm,(InterpFrm(j),j=1,PolyOrder+1)
c110  &  format(1x,f10.5,10i5)
c
        do mrk= 1, MP.num_mrk
          do axis= x, z
            do j= 1, PolyOrder+1
              XData(j)= MrkData(axis,InterpFrm(j),mrk)
            end do
            call lagrange_interp( YData,XData,PolyOrder+1,NewTime,
     &                            NewMrkData(axis,NewPnt,mrk),error )
            if ( error ) goto 2000
          end do ! axis
        end do ! mrk
        NewTime= NewTime + TimeInc
      end do ! NewPnt
      do frm= 1, NumPnts + NumPadPnts*2
        do mrk= 1, MP.num_mrk                                
          do axis= x, z
            MrkData(axis,frm,mrk)= NewMrkData(axis,frm,mrk)
          end do
        end do
      end do

c                                      
c     Update the various quantities describing the time interval of the data
c
      FrmRateRatio= NewFrmRate/rMtnFrmRate
      MP.first_frm= 1
      if ( UseCycleInfo ) then
        FirstFrm= MP.cycle_first_frm(1)
        do i= 1, MP.num_rhs
          MP.rhs_frm(i)= float(MP.rhs_frm(i)-FirstFrm+1)*
     &                        FrmRateRatio + NumPadPnts
          if ( MP.rhs_frm(i) .le. NumPadPnts ) then
            MP.rhs_frm(i)= NumPadPnts + 1
          else if ( MP.rhs_frm(i) .gt. (NumPnts + NumPadPnts) ) then
            MP.rhs_frm(i)= NumPnts + NumPadPnts
          end if
        end do  
        do i= 1, MP.num_lhs
          MP.lhs_frm(i)= float(MP.lhs_frm(i)-FirstFrm+1)*
     &                        FrmRateRatio + NumPadPnts
          if ( MP.lhs_frm(i) .le. NumPadPnts ) then
            MP.lhs_frm(i)= NumPadPnts + 1
          else if ( MP.lhs_frm(i) .gt. (NumPnts + NumPadPnts) ) then
            MP.lhs_frm(i)= NumPnts + NumPadPnts
          end if
        end do           
        do i= 1, MP.num_rto
          MP.rto_frm(i)= float(MP.rto_frm(i)-FirstFrm+1)*
     &                        FrmRateRatio + NumPadPnts
          if ( MP.rto_frm(i) .le. NumPadPnts ) then
            MP.rto_frm(i)= NumPadPnts + 1
          else if ( MP.rto_frm(i) .gt. (NumPnts + NumPadPnts) ) then
            MP.rto_frm(i)= NumPnts + NumPadPnts
          end if
        end do  
        do i= 1, MP.num_lto
          MP.lto_frm(i)= float(MP.lto_frm(i)-FirstFrm+1)*
     &                        FrmRateRatio + NumPadPnts
          if ( MP.lto_frm(i) .le. NumPadPnts ) then
            MP.lto_frm(i)= NumPadPnts + 1
          else if ( MP.lto_frm(i) .gt. (NumPnts + NumPadPnts) ) then
            MP.lto_frm(i)= NumPnts + NumPadPnts
          end if
        end do  
        if ( .not. CheckGaitPattern(OF) ) goto 2030 ! error
        MP.first_frm= FirstFrm*FrmRateRatio - NumPadPnts
      end if ! UseCycleInfo
      MP.nfrm= NumPnts + 2*NumPadPnts                       
c      MP.frm_rate= jnint(NewFrmRate)
      rMtnFrmRate= NewFrmRate
c
c     Lastly invalidate all other kinematic calculations because they now
c     are invalid since a different set of marker data is present
c
      MF.diagnostics=  .false.
      MF.VelAcc=       .false.
      
      SF.position=     .false.
      SF.vel_acc=      .false.    
      SF.energy=       .false.
      SF.FiniteScrew=  .false.
      SF.InstantScrew= .false.
      SF.cntr_mass=    .false.
      SF.momentum=     .false.
      SF.angle=        .false.
      SF.adjusted=     .false.
      SF.RefCorrect=   .false.
      SF.BCMVel=       .false.

      JF.angle_JCS=     .false.
      JF.vel_acc_JCS=   .false.
      JF.angle_EUL=     .false.
      JF.vel_acc_EUL=   .false.
      JF.load=          .false.
      JF.wrench=        .false.
      JF.power=         .false.
      JF.FiniteScrew=   .false.
      JF.InstantScrew=  .false.
      JF.legsummed=     .false.
      JF.TRU_angle=     .false.
      JF.Translation=   .false.
      JF.ScrewAngle=    .false.
      JF.MomentArm=     .false.
      JF.WrenchArm=     .false.
      JF.RefCorrectAng= .false.

      MF.Regularized= .true.
      if ( MF.NumProcOps .le. 4 ) then
        MF.NumProcOps= MF.NumProcOps + 1
        MF.ProcOp(MF.NumProcOps)= CMND.Text(1)(1:40)
      end if
      return                  

1000  write(6,1005) 
      if ( OF.ToLog ) write(3,1005)
1005  format(' Error- Marker data cannot be regularized until it has first'/
     &       '        been adjusted (MRK ADJ)')
      return

1010  write(6,1015) 
      if ( OF.ToLog ) write(3,1015)
1015  format(' Error- Marker data cannot be regularized with the gait cycle'/
     &       '        information because it has not been input yet')
      return                         
                                              
1020  write(6,1025) StartFrm,MP.nfrm+MP.first_frm-1
      if ( OF.ToLog ) write(3,1025) StartFrm,
     &                                      MP.nfrm+MP.first_frm-1
1025  format(' Error- ',i4,' is not a valid starting frame number'/
     &       '        Valid numbers are from 1 to ',i4)
      return

1030  write(6,1035) EndFrm,StartFrm,MP.nfrm+MP.first_frm-1
      if ( OF.ToLog ) write(3,1035) EndFrm,StartFrm,
     &                                      MP.nfrm+MP.first_frm-1
1035  format(' Error- ',i4,' is not a valid ending frame number'/
     &       '        Valid numbers are from ',i4,' to ',i4)
      return
                          
1040  write(6,1045) CMND.OptText(i,2),MAX_MTN
      if ( OF.ToLog ) write(3,1045) CMND.OptText(i,2),MAX_MTN
1045  format(' Error- ',a3,' is not a valid number of interpolation points'/
     &       '        Valid numbers are from 1 to ',i3)
      return
                                                    
1050  write(6,1055) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,1055) CMND.OptText(i,2)
1055  format(' Error- ',a3,' is not a valid MRK REG option'/
     &       '        Valid options are: CYCle STArt     END '/
     &       '                           PAD   POLyorder PoiNT')
      return
                                     
1060  write(6,1065) 
      if ( OF.ToLog ) write(3,1065)
1065  format(' Error- Starting frame number cannot be set if gait cycle'/
     &       '        information is to be used')
      return

1070  write(6,1075) 
      if ( OF.ToLog ) write(3,1075)
1075  format(' Error- Ending frame number cannot be set if gait cycle'/
     &       '        information is to be used')
      return

1080  write(6,1085) rMtnFrmRate
      if ( OF.ToLog ) write(3,1085) rMtnFrmRate
1085  format(' Error- The marker data frame rate is ',f10.4,' which is not'/
     &       '        valid. An internal program error has occured')
      return

2000  write(6,2005) mrk,newpnt,OldFrm
      if ( OF.ToLog ) write(3,2005) mrk,newpnt,OldFrm
2005  format(' Error- An error has occured during interpolation of marker ',i2/
     &       '        at new frame number ',i3,' (Old frame # ',i3,')')
      return

2010  write(6,2015) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,2015) CMND.OptText(i,2)             
2015  format(' Error- ',a3,' is not a valid polynomial order'/
     &       '        Valid orders are from 1 to 9')
      return

2020  write(6,2025) CMND.OptText(i,2)
      if ( OF.ToLog ) write(3,2025) CMND.OptText(i,2)             
2025  format(' Error- ',a3,' is not a valid number of start/end pad points'/
     &       '        Valid numbers are from 0 to 10')
      return

2030  write(6,2035)
      if ( OF.ToLog ) write(3,2035)
2035  format(' Error- After regularization using the original gait cycle information'/
     &       '        an error occured in updating the new modified gait cycle')
      return
      end
      
c                                        
c	MrkPosSegLCS
c
c	Copyright by Dwight Meglan 1990
c
c	Purpose:	Compute the position of the markers attched to a segment 
c           relative to the embedded local coordinate system of the segment
c           Also, calculate the relationship between the technical CS and
c           the anatomic LCS of the segment. These are stored in SP
c
c
c	Routines called:
c
c	Common Block Variables USED by routine:
c	Block Name		Variable Name
c
c
c	Common Block Variables MODIFIED by routine:
c	Block Name		Variable Name
c
      subroutine MrkPosSegLCS( CMND,MF,SF,OF )
      implicit none

      include 'a_const.inc'
      include 'a_struct.inc'
      include 'a_prm.inc'
      include 'a_names.inc'
      include 'a_mrk.inc'
                                                             
      record /CmndLine/     CMND
      record /marker_flag/  MF
      record /segment_flag/ SF
      record /option_flag/  OF

      if ( (GP.Analysis_Type .eq. CCR_ANL) .or.
     &     (GP.Analysis_Type .eq. CCL_ANL)      ) then
        call CalcCCAlignment( CMND,MF,SF,OF )
                                                   
      else if ( (GP.Analysis_Type .eq. MIR_ANL) .or.
     &          (GP.Analysis_Type .eq. MIL_ANL)      ) then
        call CalcMayoAlignment( CMND,MF,SF,OF )

      else if ( (GP.Analysis_Type .eq. F3DR_ANL) .or.
     &          (GP.Analysis_Type .eq. F3DL_ANL)      ) then
        call CalcF3DAlignment( CMND,MF,SF,OF )

      else
        write(6,100) 
        if ( OF.ToLog ) write(3,100)
100     format(' Error- Unsupported marker set'/
     &         '        Supported marker sets are: '/
     &         '           Modified Cleveland Clinic right and left leg alignment'/
     &         '           Right and left Foot3d alignment')
      end if
      return
      end
