************************************************************************
*                                                                      *
*    Data output for anybody subroutines                               *
*                                                                      *
************************************************************************


*****************
* body geometry *
*****************

c      subroutine anygeometry(len(1),len(2),len(3),len(4),len(5),len(6)
c     &                ,len(7),len(10),len(12),len(14),len(16),len(18)
c     &                ,bu(5),bu(6),bu(8),bu(9))

       subroutine anygeometry(len,bu)

      implicit double precision(a-z)
      integer j,numb
      

       parameter(numb = 21)
       parameter(j = 9)

      double precision size1,size2,size3,size4,size5,size6,size7
     &                 ,len(numb),bu(j)
c     &                 ,len(1),len(2),len(3),len(4),len(5),len(6)
c     &               ,len(7),len(10),len(12),len(14),len(16),len(18)
c     &                ,bu(5),bu(6),bu(8),bu(9)

      size1 = 0.0d0
      size2 = 0.0d0
      size3 = 0.0d0
      size4 = 0.0d0
      size5 = 0.0d0
      size6 = 0.0d0
      size7 = 0.0d0


      size1 = (len(10)+bu(5))*bu(8)
      size2 = (len(12)+bu(6))*bu(8)
      size3 = len(14)*bu(8)
      size4 = 0.16d0*bu(8)/1.759d0
      size5 = (len(1)+len(2)+len(3)+len(4)+len(5)+len(6)
     &                +len(7))*bu(8) - 0.14d0
      size6 = len(16)*bu(8)
      size7 = len(18)*bu(8)


       open(unit=10,file='AnyBody/BodyGeometry/ThighLength.any'
     &            ,status='unknown',access='append')
      write(10,*) size1
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/ShankLength.any'
     &            ,status='unknown',access='append')
      write(10,*) size2
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/FootLength.any'
     &            ,status='unknown',access='append')
      write(10,*) size3
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/PelvisWidth.any'
     &            ,status='unknown',access='append')
      write(10,*) size4
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/HeadHeight.any'
     &            ,status='unknown',access='append')
      write(10,*) 0.14
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/TrunkHeight.any'
     &            ,status='unknown',access='append')
      write(10,*) size5
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/UpperArmLength.any'
     &            ,status='unknown',access='append')
      write(10,*) size6
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/LowerArmLength.any'
     &            ,status='unknown',access='append')
      write(10,*) size7
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/BodyHeight.any'
     &            ,status='unknown',access='append')
      write(10,*) bu(8)
         close(10)

       open(unit=10,file='AnyBody/BodyGeometry/BodyMass.any'
     &            ,status='unknown',access='append')
      write(10,*) bu(9)
         close(10)

      return
      end





****************
* joint motion *
****************
      subroutine outanyjoint(numb,dummy,jcount,pi,ttt)

      implicit double precision(a-z)
      integer numb,dummy,jcount,jseg(jcount),jax(jcount)
     &            ,jjj,jjk,jfr,jfr2,jfr3,frame,conv(jcount)
     &            ,jjl,jjm,jjn,jjo,jjp,jjq
      
      double precision pi,ttt
     &                ,tFrame
     &                ,deg(jcount,dummy),dega1(numb,dummy)
     &                ,dega2(numb,dummy),dega3(numb,dummy)
     &                ,ghr1x(dummy),ghr1y(dummy),ghr1z(dummy)
     &                ,ghr2x(dummy),ghr2y(dummy),ghr2z(dummy)
     &                ,ghl1x(dummy),ghl1y(dummy),ghl1z(dummy)
     &                ,ghl2x(dummy),ghl2y(dummy),ghl2z(dummy)
     &                ,hpr1x(dummy),hpr1y(dummy),hpr1z(dummy)
     &                ,hpr2x(dummy),hpr2y(dummy),hpr2z(dummy)
     &                ,hpl1x(dummy),hpl1y(dummy),hpl1z(dummy)
     &                ,hpl2x(dummy),hpl2y(dummy),hpl2z(dummy)


      do 517 jjm=1,jcount
      do 516 jjl=1,dummy
      deg(jjm,jjl)=0.0d0
  516 continue
  517 continue

      do 519 jjn=1,numb
      do 518 jjo=1,dummy
      dega1(jjn,jjo) = 0.0d0
      dega2(jjn,jjo) = 0.0d0
      dega3(jjn,jjo) = 0.0d0
  518 continue
  519 continue

      do 520 jjp=1,dummy
      ghr1x(jjp) = 0.0d0
      ghr1y(jjp) = 0.0d0
      ghr1z(jjp) = 1.0d0

      ghr2x(jjp) = 1.0d0
      ghr2y(jjp) = 0.0d0
      ghr2z(jjp) = 1.0d0

      ghl1x(jjp) = 0.0d0
      ghl1y(jjp) = 0.0d0
      ghl1z(jjp) = 1.0d0

      ghl2x(jjp) = 1.0d0
      ghl2y(jjp) = 0.0d0
      ghl2z(jjp) = 1.0d0

      hpr1x(jjp) = 0.0d0
      hpr1y(jjp) = 0.0d0
      hpr1z(jjp) = -1.0d0

      hpr2x(jjp) = 1.0d0
      hpr2y(jjp) = 0.0d0
      hpr2z(jjp) = -1.0d0

      hpl1x(jjp) = 0.0d0
      hpl1y(jjp) = 0.0d0
      hpl1z(jjp) = -1.0d0

      hpl2x(jjp) = 1.0d0
      hpl2y(jjp) = 0.0d0
      hpl2z(jjp) = -1.0d0
  520 continue

***** Set initial condition of SWUM ******************
      do 521 jjq=1,dummy
      dega2(7,jjq) = dega2(7,jjq) + 20.0d0
      dega1(5,jjq) = dega1(5,jjq) - 20.0d0
      dega3(5,jjq) = dega3(5,jjq) - 20.0d0
      dega3(18,jjq) = dega3(18,jjq) - 30.0d0
      dega3(19,jjq) = dega3(19,jjq) - 30.0d0
      dega2(14,jjq) = dega2(14,jjq) - 90.0d0
      dega2(15,jjq) = dega2(15,jjq) - 90.0d0
  521 continue

 
      open(unit=24,file='joint_motion.dat',status='unknown')
      read(24,*) frame

c      do 510 jjj=1,jcount+1
       do 511 jjj=1,10000

        read(24,*) jseg(jjj),jax(jjj)

c       if (jjj.eq.jcount+1) then
c          goto 225
c       end if
      if (jseg(jjj).eq.0) then
         goto 225
      end if

      do 512 jfr=1,frame
         read(24,*) deg(jjj,jfr)
  512 continue
c  510 continue
  511 continue
  225 continue




***************recognition***********************************************

      do 513 jjk=1,jcount
       read(24,*) conv(jjk)

      if (conv(jjk).eq.1) then

      do 514 jfr2=1,frame


       if (jseg(jjk).eq.1) then
********************
* (1) Lower waist  *
********************
       if (jax(jjk).eq.1) then
         dega1(1,jfr2) = dega1(1,jfr2) + deg(jjk,jfr2)
       
       else if (jax(jjk).eq.2) then
         dega2(1,jfr2) = dega2(1,jfr2) - deg(jjk,jfr2) 

       else if (jax(jjk).eq.3) then
         dega3(1,jfr2) = dega3(1,jfr2) + deg(jjk,jfr2) 

        end if


      else if (jseg(jjk).eq.4) then
********************
* (4) upper breast *
********************
       if (jax(jjk).eq.4) then
         dega1(4,jfr2) = dega1(4,jfr2) - deg(jjk,jfr2)
       
       else if (jax(jjk).eq.7) then
         dega2(4,jfr2) = dega2(4,jfr2) + deg(jjk,jfr2) 

       else if (jax(jjk).eq.6) then
         dega1(5,jfr2) = dega1(5,jfr2) + deg(jjk,jfr2) 

       else if (jax(jjk).eq.9) then
         dega3(5,jfr2) = dega3(5,jfr2) - deg(jjk,jfr2)  
        end if


      else if (jseg(jjk).eq.3) then
********************
* (3) lower breast *
********************
       if (jax(jjk).eq.1) then
         dega1(3,jfr2) = dega1(3,jfr2) + deg(jjk,jfr2)
       
       else if (jax(jjk).eq.2) then
         dega2(3,jfr2) = dega2(3,jfr2) - deg(jjk,jfr2) 

       else if (jax(jjk).eq.3) then
         dega3(3,jfr2) = dega3(3,jfr2) + deg(jjk,jfr2) 

        end if


      else if (jseg(jjk).eq.7) then
********************
* (7)  head        *
********************
       if (jax(jjk).eq.2) then
         dega2(7,jfr2) = dega2(7,jfr2) - deg(jjk,jfr2)
        end if 


      else if (jseg(jjk).eq.18) then
************************
* (18) right forearm *
************************
        if (jax(jjk).eq.3) then
         dega3(18,jfr2) = dega3(18,jfr2) - deg(jjk,jfr2) 

        else if (jax(jjk).eq.2) then
         dega2(18,jfr2) = dega2(18,jfr2) - deg(jjk,jfr2) 
        end if 


      else if (jseg(jjk).eq.19) then
***********************
* (19) Left forearm *
***********************
        if (jax(jjk).eq.3) then
         dega3(19,jfr2) = dega3(19,jfr2) + deg(jjk,jfr2) 

        else if (jax(jjk).eq.2) then
         dega2(19,jfr2) = dega2(19,jfr2) - deg(jjk,jfr2) 
        end if


      else if (jseg(jjk).eq.20) then
********************
* (20) Right wrist *
********************
        if (jax(jjk).eq.2) then
         dega2(20,jfr2) = dega2(20,jfr2) + deg(jjk,jfr2) 

        else if (jax(jjk).eq.1) then
         dega1(20,jfr2) = dega1(20,jfr2) + deg(jjk,jfr2) 
        end if 


       else if (jseg(jjk).eq.21) then
********************
* (21) Left wrist *
********************
        if (jax(jjk).eq.2) then
         dega2(21,jfr2) = dega2(21,jfr2) + deg(jjk,jfr2) 

        else if (jax(jjk).eq.1) then
         dega1(21,jfr2) = dega1(21,jfr2) - deg(jjk,jfr2) 
        end if 


       else if (jseg(jjk).eq.12) then
************************
* (12) right shank     *
************************
        if (jax(jjk).eq.2) then
         dega2(12,jfr2) = dega2(12,jfr2) + deg(jjk,jfr2) 

        else if (jax(jjk).eq.3) then
         dega3(12,jfr2) = dega3(12,jfr2) + deg(jjk,jfr2) 
        end if 


      else if (jseg(jjk).eq.13) then
************************
* (13) left shank      *
************************
       if (jax(jjk).eq.2) then
         dega2(13,jfr2) = dega2(13,jfr2) + deg(jjk,jfr2) 

       else if (jax(jjk).eq.3) then
         dega3(13,jfr2) = dega3(13,jfr2) - deg(jjk,jfr2) 
        end if 


      else if (jseg(jjk).eq.14) then
************************
* (14) right ankle     *
************************
        if (jax(jjk).eq.2) then
         dega2(14,jfr2) = dega2(14,jfr2) - deg(jjk,jfr2) 

        else if (jax(jjk).eq.1) then
         dega1(14,jfr2) = dega1(14,jfr2) - deg(jjk,jfr2) 
        end if 

      else if (jseg(jjk).eq.15) then
************************
* (15) left ankle      *
************************
        if (jax(jjk).eq.2) then
         dega2(15,jfr2) = dega2(15,jfr2) - deg(jjk,jfr2)

        else if (jax(jjk).eq.1) then
         dega1(15,jfr2) = dega1(15,jfr2) + deg(jjk,jfr2) 
        end if 


      else if (jseg(jjk).eq.16) then
************************
* (16) right upper arm *
************************
        deg(jjk,jfr2) = deg(jjk,jfr2)*pi/180.0d0

        if (jax(jjk).eq.1) then
         call vecrot(ghr1x(jfr2),ghr1y(jfr2),ghr1z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)
         call vecrot(ghr2x(jfr2),ghr2y(jfr2),ghr2z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)

      else if (jax(jjk).eq.2) then
         call vecrot(ghr1x(jfr2),ghr1y(jfr2),ghr1z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)
         call vecrot(ghr2x(jfr2),ghr2y(jfr2),ghr2z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)

      else if (jax(jjk).eq.3) then
         call vecrot(ghr1x(jfr2),ghr1y(jfr2),ghr1z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
         call vecrot(ghr2x(jfr2),ghr2y(jfr2),ghr2z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
      end if


      else if (jseg(jjk).eq.17) then
************************
* (17) left upper arm  *
************************
        deg(jjk,jfr2) = deg(jjk,jfr2)*pi/180.0d0

        if (jax(jjk).eq.1) then
         call vecrot(ghl1x(jfr2),ghl1y(jfr2),ghl1z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)
         call vecrot(ghl2x(jfr2),ghl2y(jfr2),ghl2z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)

      else if (jax(jjk).eq.2) then
         call vecrot(ghl1x(jfr2),ghl1y(jfr2),ghl1z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)
         call vecrot(ghl2x(jfr2),ghl2y(jfr2),ghl2z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)

      else if (jax(jjk).eq.3) then
         call vecrot(ghl1x(jfr2),ghl1y(jfr2),ghl1z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
         call vecrot(ghl2x(jfr2),ghl2y(jfr2),ghl2z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
      end if


      else if (jseg(jjk).eq.10) then
************************
* (10) right thigh     *
************************
        deg(jjk,jfr2) = deg(jjk,jfr2)*pi/180.0d0

        if (jax(jjk).eq.1) then
         call vecrot(hpr1x(jfr2),hpr1y(jfr2),hpr1z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)
         call vecrot(hpr2x(jfr2),hpr2y(jfr2),hpr2z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)

      else if (jax(jjk).eq.2) then
         call vecrot(hpr1x(jfr2),hpr1y(jfr2),hpr1z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)
         call vecrot(hpr2x(jfr2),hpr2y(jfr2),hpr2z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)

      else if (jax(jjk).eq.3) then
         call vecrot(hpr1x(jfr2),hpr1y(jfr2),hpr1z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
         call vecrot(hpr2x(jfr2),hpr2y(jfr2),hpr2z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
      end if


      else if (jseg(jjk).eq.11) then
************************
* (11) left thigh      *
************************
        deg(jjk,jfr2) = deg(jjk,jfr2)*pi/180.0d0

        if (jax(jjk).eq.1) then
         call vecrot(hpl1x(jfr2),hpl1y(jfr2),hpl1z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)
         call vecrot(hpl2x(jfr2),hpl2y(jfr2),hpl2z(jfr2)
     &                               ,deg(jjk,jfr2),0.0d0,0.0d0)

      else if (jax(jjk).eq.2) then
         call vecrot(hpl1x(jfr2),hpl1y(jfr2),hpl1z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)
         call vecrot(hpl2x(jfr2),hpl2y(jfr2),hpl2z(jfr2)
     &                               ,0.0d0,deg(jjk,jfr2),0.0d0)

      else if (jax(jjk).eq.3) then
         call vecrot(hpl1x(jfr2),hpl1y(jfr2),hpl1z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
         call vecrot(hpl2x(jfr2),hpl2y(jfr2),hpl2z(jfr2)
     &                               ,0.0d0,0.0d0,deg(jjk,jfr2))
      end if


      end if

  514 continue
       end if

  513 continue

      close(24)



********output************************************************

c 9002 format(1p6e17.8)
      do 515 jfr3=1,frame

************************
* (4)  upper breast    *
************************
       open(unit=10,file='AnyBody/JointMotion/SCElevationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(4,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(4,1)
           end if
*------------------------------------
         close(10)

       open(unit=10,file='AnyBody/JointMotion/SCElevationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(4,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(4,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/SCAxialRotationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(5,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(5,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/SCAxialRotationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(5,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(5,1)
           end if
*------------------------------------
         close(10)


      open(unit=10,file='AnyBody/JointMotion/SCProtractionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) '-23.0,'
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) -23.0
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/SCProtractionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) '-23.0,'
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) -23.0
           end if
*------------------------------------
         close(10)


************************
* (3)  lower breast    *
************************
       open(unit=10,file='AnyBody/JointMotion/PTExtension.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(3,jfr3)-dega2(1,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(3,1)-dega2(1,1)
           end if
*------------------------------------
         close(10)

       open(unit=10,file='AnyBody/JointMotion/PTLateralBending.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(3,jfr3)-dega1(1,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(3,1)-dega1(1,1)
           end if
*------------------------------------
         close(10)

       open(unit=10,file='AnyBody/JointMotion/PTRotation.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(3,jfr3)-dega3(1,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(3,1)-dega3(1,1)
           end if
*------------------------------------
         close(10)


************************
* (7)  head             *
************************
       open(unit=10,file='AnyBody/JointMotion/NeckExtension.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(7,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(7,1)
           end if
*------------------------------------
         close(10)

************************
* (18) right forearm   *
************************
       open(unit=10,file='AnyBody/JointMotion/ElbowPronationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(18,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(18,1)
           end if
*------------------------------------
        close(10)

       open(unit=10,file='AnyBody/JointMotion/ElbowFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(18,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(18,1)
           end if
*------------------------------------
         close(10)

***********************
* (19) Left forearm   *
***********************
      open(unit=10,file='AnyBody/JointMotion/ElbowPronationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(19,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(19,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/ElbowFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(19,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(19,1)
           end if
*------------------------------------
         close(10)

********************
* (20) Right wrist *
********************
      open(unit=10,file='AnyBody/JointMotion/WristFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(20,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(20,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/WristAbductionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(20,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(20,1)
           end if
*------------------------------------
         close(10)

********************
* (21) Left wrist *
********************

      open(unit=10,file='AnyBody/JointMotion/WristFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(21,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(21,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/WristAbductionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(21,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(21,1)
           end if
*------------------------------------
         close(10)

************************
* (12) right shank     *
************************

      open(unit=10,file='AnyBody/JointMotion/KneeFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(12,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(12,1)
           end if
*------------------------------------
         close(10)

      open(unit=10
     &       ,file='AnyBody/JointMotion/KneeExternalRotationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(12,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(12,1)
           end if
*------------------------------------
         close(10)


      open(unit=10,file='AnyBody/JointMotion/KneeAbductionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(12,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(12,1)
           end if
*------------------------------------
         close(10)

     


************************
* (13) left shank      *
************************

      open(unit=10,file='AnyBody/JointMotion/KneeFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(13,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(13,1)
           end if
*------------------------------------
         close(10)

      open(unit=10
     &     ,file='AnyBody/JointMotion/KneeExternalRotationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(13,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega3(13,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/KneeAbductionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(13,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(13,1)
           end if
*------------------------------------
         close(10)


************************
* (14) right ankle     *
************************

      open(unit=10
     &       ,file='AnyBody/JointMotion/AnklePlantarFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(14,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(14,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/AnkleEversionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(14,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(14,1)
           end if
*------------------------------------
         close(10)

************************
* (15) left ankle      *
************************

      open(unit=10
     &      ,file='AnyBody/JointMotion/AnklePlantarFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(15,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega2(15,1)
           end if
*------------------------------------
         close(10)

      open(unit=10,file='AnyBody/JointMotion/AnkleEversionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(15,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then
            write(10,*) dega1(15,1)
           end if
*------------------------------------
         close(10)



************************
* (16) right upper arm *
************************


*----Flexion------*


      if ((ghr1x(jfr3).eq.0.0d0).and.(ghr1z(jfr3).ge.0.0d0)) then
            dega2(16,jfr3) = 0.0d0
      else if ((ghr1x(jfr3).eq.0.0d0).and.(ghr1z(jfr3).lt.0.0d0)) then
            dega2(16,jfr3) = pi
      else
       dega2(16,jfr3) = dacos(ghr1z(jfr3)
     &                   /dsqrt(ghr1z(jfr3)**2+ghr1x(jfr3)**2))
     &                  + (1-ghr1x(jfr3)/dabs(ghr1x(jfr3)))
     &                   * (pi-dacos(ghr1z(jfr3)
     &                    /dsqrt(ghr1z(jfr3)**2+ghr1x(jfr3)**2)))
      end if


*----Abduction-----*
         call vecrot(ghr1x(jfr3),ghr1y(jfr3),ghr1z(jfr3)
     &                             ,0.0d0,-dega2(16,jfr3),0.0d0)
         call vecrot(ghr2x(jfr3),ghr2y(jfr3),ghr2z(jfr3)
     &                             ,0.0d0,-dega2(16,jfr3),0.0d0)

      if ((ghr1z(jfr3).eq.0.0d0).and.(ghr1y(jfr3).ge.0.0d0)) then
            dega1(16,jfr3) = 0.0d0 - pi/2
      else if ((ghr1z(jfr3).eq.0.0d0).and.(ghr1y(jfr3).lt.0.0d0)) then
            dega1(16,jfr3) = pi - pi/2
      else
       dega1(16,jfr3) = dacos(ghr1y(jfr3)
     &                   /dsqrt(ghr1y(jfr3)**2+ghr1z(jfr3)**2))
     &                  + (1-ghr1z(jfr3)/dabs(ghr1z(jfr3)))
     &                   * (pi-dacos(ghr1y(jfr3)
     &                    /dsqrt(ghr1y(jfr3)**2+ghr1z(jfr3)**2)))
     &                  - pi/2
      end if

*---External Rotation---*
         call vecrot(ghr1x(jfr3),ghr1y(jfr3),ghr1z(jfr3)
     &                             ,-dega1(16,jfr3),0.0d0,0.0d0)
         call vecrot(ghr2x(jfr3),ghr2y(jfr3),ghr2z(jfr3)
     &                             ,-dega1(16,jfr3),0.0d0,0.0d0)

      if ((ghr2y(jfr3).eq.0.0d0).and.(ghr2x(jfr3).ge.0.0d0)) then
            dega3(16,jfr3) = 0.0d0
      else if ((ghr2y(jfr3).eq.0.0d0).and.(ghr2x(jfr3).lt.0.0d0)) then
            dega3(16,jfr3) = pi
      else
       dega3(16,jfr3) = dacos(ghr2x(jfr3)
     &                   /dsqrt(ghr2x(jfr3)**2+ghr2y(jfr3)**2))
     &                  + (1-ghr2y(jfr3)/dabs(ghr2y(jfr3)))
     &                   * (pi-dacos(ghr2x(jfr3)
     &                    /dsqrt(ghr2x(jfr3)**2+ghr2y(jfr3)**2)))
      end if

******Convert for Anybody****************************************

      dega2(16,jfr3) = -dega2(16,jfr3)*180.0d0/pi-dega2(3,jfr3)

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega2(16,jfr3)-dega2(16,jfr3-1)).gt.180.0d0)) then
            dega2(16,jfr3) = dega2(16,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega2(16,jfr3)-dega2(16,jfr3-1)).lt.-180.0d0)) then
            dega2(16,jfr3) = dega2(16,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega1(16,jfr3) = -dega1(16,jfr3)*180.0d0/pi-dega1(3,jfr3)

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega1(16,jfr3)-dega1(16,jfr3-1)).gt.180.0d0)) then
            dega1(16,jfr3) = dega1(16,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega1(16,jfr3)-dega1(16,jfr3-1)).lt.-180.0d0)) then
            dega1(16,jfr3) = dega1(16,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega3(16,jfr3) = -dega3(16,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega3(16,jfr3)-dega3(16,jfr3-1)).gt.180.0d0)) then
            dega3(16,jfr3) = dega3(16,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega3(16,jfr3)-dega3(16,jfr3-1)).lt.-180.0d0)) then
            dega3(16,jfr3) = dega3(16,jfr3) + 360.0d0
      end if 
*-------------------------------------

******Output to file****************************************
      open(unit=10,file='AnyBody/JointMotion/GHFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(16,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega2(16,1)-dega2(16,jfr3)).gt.180.0d0)) then
             dega2(16,1) = dega2(16,1) - 360.0d0
            else if (((dega2(16,1)-dega2(16,jfr3)).lt.-180.0d0)) then
             dega2(16,1) = dega2(16,1) + 360.0d0
            end if 

            write(10,*) dega2(16,1)
           end if
*------------------------------------
        close(10)

      open(unit=10,file='AnyBody/JointMotion/GHAbductionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(16,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega1(16,1)-dega1(16,jfr3)).gt.180.0d0)) then
             dega1(16,1) = dega1(16,1) - 360.0d0
            else if (((dega1(16,1)-dega1(16,jfr3)).lt.-180.0d0)) then
             dega1(16,1) = dega1(16,1) + 360.0d0
            end if

            write(10,*) dega1(16,1)
           end if
*------------------------------------
        close(10)

      open(unit=10
     &        ,file='AnyBody/JointMotion/GHExternalRotationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(16,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega3(16,1)-dega3(16,jfr3)).gt.180.0d0)) then
             dega3(16,1) = dega3(16,1) - 360.0d0
            else if (((dega3(16,1)-dega3(16,jfr3)).lt.-180.0d0)) then
             dega3(16,1) = dega3(16,1) + 360.0d0
            end if

            write(10,*) dega3(16,1)
           end if
*------------------------------------
        close(10)


************************
* (17) left upper arm  *
************************


*----Flexion------*

      if ((ghl1x(jfr3).eq.0.0d0).and.(ghl1z(jfr3).ge.0.0d0)) then
            dega2(17,jfr3) = 0.0d0
      else if ((ghl1x(jfr3).eq.0.0d0).and.(ghl1z(jfr3).lt.0.0d0)) then
            dega2(17,jfr3) = pi
      else
       dega2(17,jfr3) = dacos(ghl1z(jfr3)
     &                   /dsqrt(ghl1z(jfr3)**2+ghl1x(jfr3)**2))
     &                  + (1-ghl1x(jfr3)/dabs(ghl1x(jfr3)))
     &                   * (pi-dacos(ghl1z(jfr3)
     &                    /dsqrt(ghl1z(jfr3)**2+ghl1x(jfr3)**2)))
      end if

*----Abduction-----*
         call vecrot(ghl1x(jfr3),ghl1y(jfr3),ghl1z(jfr3)
     &                             ,0.0d0,-dega2(17,jfr3),0.0d0)
         call vecrot(ghl2x(jfr3),ghl2y(jfr3),ghl2z(jfr3)
     &                             ,0.0d0,-dega2(17,jfr3),0.0d0)

      if ((ghl1z(jfr3).eq.0.0d0).and.(ghl1y(jfr3).ge.0.0d0)) then
            dega1(17,jfr3) = 0.0d0 - pi/2
      else if ((ghl1z(jfr3).eq.0.0d0).and.(ghl1y(jfr3).lt.0.0d0)) then
            dega1(17,jfr3) = pi - pi/2
      else
       dega1(17,jfr3) = dacos(ghl1y(jfr3)
     &                   /dsqrt(ghl1y(jfr3)**2+ghl1z(jfr3)**2))
     &                  + (1-ghl1z(jfr3)/dabs(ghl1z(jfr3)))
     &                   * (pi-dacos(ghl1y(jfr3)
     &                    /dsqrt(ghl1y(jfr3)**2+ghl1z(jfr3)**2)))
     &                  - pi/2
      end if

*---External Rotation---*
         call vecrot(ghl1x(jfr3),ghl1y(jfr3),ghl1z(jfr3)
     &                             ,-dega1(17,jfr3),0.0d0,0.0d0)
         call vecrot(ghl2x(jfr3),ghl2y(jfr3),ghl2z(jfr3)
     &                             ,-dega1(17,jfr3),0.0d0,0.0d0)

      if ((ghl2y(jfr3).eq.0.0d0).and.(ghl2x(jfr3).ge.0.0d0)) then
            dega3(17,jfr3) = 0.0d0
      else if ((ghl2y(jfr3).eq.0.0d0).and.(ghl2x(jfr3).lt.0.0d0)) then
            dega3(17,jfr3) = pi
      else
       dega3(17,jfr3) = dacos(ghl2x(jfr3)
     &                   /dsqrt(ghl2x(jfr3)**2+ghl2y(jfr3)**2))
     &                  + (1-ghl2y(jfr3)/dabs(ghl2y(jfr3)))
     &                   * (pi-dacos(ghl2x(jfr3)
     &                    /dsqrt(ghl2x(jfr3)**2+ghl2y(jfr3)**2)))
      end if

******Convert for Anybody****************************************

      dega2(17,jfr3) = -dega2(17,jfr3)*180.0d0/pi-dega2(3,jfr3)

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega2(17,jfr3)-dega2(17,jfr3-1)).gt.180.0d0)) then
            dega2(17,jfr3) = dega2(17,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega2(17,jfr3)-dega2(17,jfr3-1)).lt.-180.0d0)) then
            dega2(17,jfr3) = dega2(17,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega1(17,jfr3) = dega1(17,jfr3)*180.0d0/pi+dega1(3,jfr3)

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega1(17,jfr3)-dega1(17,jfr3-1)).gt.180.0d0)) then
            dega1(17,jfr3) = dega1(17,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega1(17,jfr3)-dega1(17,jfr3-1)).lt.-180.0d0)) then
            dega1(17,jfr3) = dega1(17,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega3(17,jfr3) = dega3(17,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega3(17,jfr3)-dega3(17,jfr3-1)).gt.180.0d0)) then
            dega3(17,jfr3) = dega3(17,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega3(17,jfr3)-dega3(17,jfr3-1)).lt.-180.0d0)) then
            dega3(17,jfr3) = dega3(17,jfr3) + 360.0d0
      end if 
*-------------------------------------

******Output to file****************************************
      open(unit=10,file='AnyBody/JointMotion/GHFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(17,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega2(17,1)-dega2(17,jfr3)).gt.180.0d0)) then
             dega2(17,1) = dega2(17,1) - 360.0d0
            else if (((dega2(17,1)-dega2(17,jfr3)).lt.-180.0d0)) then
             dega2(17,1) = dega2(17,1) + 360.0d0
            end if 

            write(10,*) dega2(17,1)
           end if
*------------------------------------
        close(10)

      open(unit=10,file='AnyBody/JointMotion/GHAbductionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(17,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega1(17,1)-dega1(17,jfr3)).gt.180.0d0)) then
             dega1(17,1) = dega1(17,1) - 360.0d0
            else if (((dega1(17,1)-dega1(17,jfr3)).lt.-180.0d0)) then
             dega1(17,1) = dega1(17,1) + 360.0d0
            end if

            write(10,*) dega1(17,1)
           end if
*------------------------------------
        close(10)

      open(unit=10
     &     ,file='AnyBody/JointMotion/GHExternalRotationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(17,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega3(17,1)-dega3(17,jfr3)).gt.180.0d0)) then
             dega3(17,1) = dega3(17,1) - 360.0d0
            else if (((dega3(17,1)-dega3(17,jfr3)).lt.-180.0d0)) then
             dega3(17,1) = dega3(17,1) + 360.0d0
            end if

            write(10,*) dega3(17,1)
           end if
*------------------------------------
        close(10)


************************
* (10) right thigh     *
************************


*----Flexion------*

      if ((hpr1x(jfr3).eq.0.0d0).and.(hpr1z(jfr3).ge.0.0d0)) then
            dega2(10,jfr3) = 0.0d0 - pi
      else if ((hpr1x(jfr3).eq.0.0d0).and.(hpr1z(jfr3).lt.0.0d0)) then
            dega2(10,jfr3) = pi - pi
      else
       dega2(10,jfr3) = dacos(hpr1z(jfr3)
     &                   /dsqrt(hpr1z(jfr3)**2+hpr1x(jfr3)**2))
     &                  + (1-hpr1x(jfr3)/dabs(hpr1x(jfr3)))
     &                   * (pi-dacos(hpr1z(jfr3)
     &                    /dsqrt(hpr1z(jfr3)**2+hpr1x(jfr3)**2)))
     &                  - pi
      end if

*----Abduction-----*
         call vecrot(hpr1x(jfr3),hpr1y(jfr3),hpr1z(jfr3)
     &                             ,0.0d0,-dega2(10,jfr3),0.0d0)
         call vecrot(hpr2x(jfr3),hpr2y(jfr3),hpr2z(jfr3)
     &                             ,0.0d0,-dega2(10,jfr3),0.0d0)

      if ((hpr1z(jfr3).eq.0.0d0).and.(hpr1y(jfr3).ge.0.0d0)) then
            dega1(10,jfr3) = 0.0d0 - pi/2 - pi
      else if ((hpr1z(jfr3).eq.0.0d0).and.(hpr1y(jfr3).lt.0.0d0)) then
            dega1(10,jfr3) = pi - pi/2 - pi
      else
       dega1(10,jfr3) = dacos(hpr1y(jfr3)
     &                   /dsqrt(hpr1y(jfr3)**2+hpr1z(jfr3)**2))
     &                  + (1-hpr1z(jfr3)/dabs(hpr1z(jfr3)))
     &                   * (pi-dacos(hpr1y(jfr3)
     &                    /dsqrt(hpr1y(jfr3)**2+hpr1z(jfr3)**2)))
     &                  - pi/2 - pi
      end if

*---External Rotation---*
         call vecrot(hpr1x(jfr3),hpr1y(jfr3),hpr1z(jfr3)
     &                             ,-dega1(10,jfr3),0.0d0,0.0d0)
         call vecrot(hpr2x(jfr3),hpr2y(jfr3),hpr2z(jfr3)
     &                             ,-dega1(10,jfr3),0.0d0,0.0d0)

      if ((hpr2y(jfr3).eq.0.0d0).and.(hpr2x(jfr3).ge.0.0d0)) then
            dega3(10,jfr3) = 0.0d0
      else if ((hpr2y(jfr3).eq.0.0d0).and.(hpr2x(jfr3).lt.0.0d0)) then
            dega3(10,jfr3) = pi
      else
       dega3(10,jfr3) = dacos(hpr2x(jfr3)
     &                   /dsqrt(hpr2x(jfr3)**2+hpr2y(jfr3)**2))
     &                  + (1-hpr2y(jfr3)/dabs(hpr2y(jfr3)))
     &                   * (pi-dacos(hpr2x(jfr3)
     &                    /dsqrt(hpr2x(jfr3)**2+hpr2y(jfr3)**2)))
      end if

******Convert for Anybody****************************************

      dega2(10,jfr3) = -dega2(10,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega2(10,jfr3)-dega2(10,jfr3-1)).gt.180.0d0)) then
            dega2(10,jfr3) = dega2(10,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega2(10,jfr3)-dega2(10,jfr3-1)).lt.-180.0d0)) then
            dega2(10,jfr3) = dega2(10,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega1(10,jfr3) = -dega1(10,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega1(10,jfr3)-dega1(10,jfr3-1)).gt.180.0d0)) then
            dega1(10,jfr3) = dega1(10,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega1(10,jfr3)-dega1(10,jfr3-1)).lt.-180.0d0)) then
            dega1(10,jfr3) = dega1(10,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega3(10,jfr3) = -dega3(10,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega3(10,jfr3)-dega3(10,jfr3-1)).gt.180.0d0)) then
            dega3(10,jfr3) = dega3(10,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega3(10,jfr3)-dega3(10,jfr3-1)).lt.-180.0d0)) then
            dega3(10,jfr3) = dega3(10,jfr3) + 360.0d0
      end if 
*-------------------------------------

******Output to file****************************************
      open(unit=10,file='AnyBody/JointMotion/HipFlexionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(10,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega2(10,1)-dega2(10,jfr3)).gt.180.0d0)) then
             dega2(10,1) = dega2(10,1) - 360.0d0
            else if (((dega2(10,1)-dega2(10,jfr3)).lt.-180.0d0)) then
             dega2(10,1) = dega2(10,1) + 360.0d0
            end if 

            write(10,*) dega2(10,1)
           end if
*------------------------------------
        close(10)

      open(unit=10,file='AnyBody/JointMotion/HipAbductionRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(10,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega1(10,1)-dega1(10,jfr3)).gt.180.0d0)) then
             dega1(10,1) = dega1(10,1) - 360.0d0
            else if (((dega1(10,1)-dega1(10,jfr3)).lt.-180.0d0)) then
             dega1(10,1) = dega1(10,1) + 360.0d0
            end if

            write(10,*) dega1(10,1)
           end if
*------------------------------------
        close(10)

      open(unit=10
     &    ,file='AnyBody/JointMotion/HipExternalRotationRight.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(10,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega3(10,1)-dega3(10,jfr3)).gt.180.0d0)) then
             dega3(10,1) = dega3(10,1) - 360.0d0
            else if (((dega3(10,1)-dega3(10,jfr3)).lt.-180.0d0)) then
             dega3(10,1) = dega3(10,1) + 360.0d0
            end if

            write(10,*) dega3(10,1)
           end if
*------------------------------------
        close(10)


************************
* (11) left thigh      *
************************


*----Flexion------*

      if ((hpl1x(jfr3).eq.0.0d0).and.(hpl1z(jfr3).ge.0.0d0)) then
            dega2(11,jfr3) = 0.0d0 - pi
      else if ((hpl1x(jfr3).eq.0.0d0).and.(hpl1z(jfr3).lt.0.0d0)) then
            dega2(11,jfr3) = pi - pi
      else
       dega2(11,jfr3) = dacos(hpl1z(jfr3)
     &                   /dsqrt(hpl1z(jfr3)**2+hpl1x(jfr3)**2))
     &                  + (1-hpl1x(jfr3)/dabs(hpl1x(jfr3)))
     &                   * (pi-dacos(hpl1z(jfr3)
     &                    /dsqrt(hpl1z(jfr3)**2+hpl1x(jfr3)**2)))
     &                  - pi
      end if

*----Abduction-----*
         call vecrot(hpl1x(jfr3),hpl1y(jfr3),hpl1z(jfr3)
     &                             ,0.0d0,-dega2(11,jfr3),0.0d0)
         call vecrot(hpl2x(jfr3),hpl2y(jfr3),hpl2z(jfr3)
     &                             ,0.0d0,-dega2(11,jfr3),0.0d0)

      if ((hpl1z(jfr3).eq.0.0d0).and.(hpl1y(jfr3).ge.0.0d0)) then
            dega1(11,jfr3) = 0.0d0 - pi/2
      else if ((hpl1z(jfr3).eq.0.0d0).and.(hpl1y(jfr3).lt.0.0d0)) then
            dega1(11,jfr3) = pi - pi/2
      else
       dega1(11,jfr3) = dacos(hpl1y(jfr3)
     &                   /dsqrt(hpl1y(jfr3)**2+hpl1z(jfr3)**2))
     &                  + (1-hpl1z(jfr3)/dabs(hpl1z(jfr3)))
     &                   * (pi-dacos(hpl1y(jfr3)
     &                    /dsqrt(hpl1y(jfr3)**2+hpl1z(jfr3)**2)))
     &                  - pi/2 - pi
      end if

*---External Rotation---*
         call vecrot(hpl1x(jfr3),hpl1y(jfr3),hpl1z(jfr3)
     &                             ,-dega1(11,jfr3),0.0d0,0.0d0)
         call vecrot(hpl2x(jfr3),hpl2y(jfr3),hpl2z(jfr3)
     &                             ,-dega1(11,jfr3),0.0d0,0.0d0)

      if ((hpl2y(jfr3).eq.0.0d0).and.(hpl2x(jfr3).ge.0.0d0)) then
            dega3(11,jfr3) = 0.0d0
      else if ((hpl2y(jfr3).eq.0.0d0).and.(hpl2x(jfr3).lt.0.0d0)) then
            dega3(11,jfr3) = pi
      else
       dega3(11,jfr3) = dacos(hpl2x(jfr3)
     &                   /dsqrt(hpl2x(jfr3)**2+hpl2y(jfr3)**2))
     &                  + (1-hpl2y(jfr3)/dabs(hpl2y(jfr3)))
     &                   * (pi-dacos(hpl2x(jfr3)
     &                    /dsqrt(hpl2x(jfr3)**2+hpl2y(jfr3)**2)))
      end if

******Convert for Anybody****************************************

      dega2(11,jfr3) = -dega2(11,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega2(11,jfr3)-dega2(11,jfr3-1)).gt.180.0d0)) then
            dega2(11,jfr3) = dega2(11,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega2(11,jfr3)-dega2(11,jfr3-1)).lt.-180.0d0)) then
            dega2(11,jfr3) = dega2(11,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega1(11,jfr3) = dega1(11,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega1(11,jfr3)-dega1(11,jfr3-1)).gt.180.0d0)) then
            dega1(11,jfr3) = dega1(11,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega1(11,jfr3)-dega1(11,jfr3-1)).lt.-180.0d0)) then
            dega1(11,jfr3) = dega1(11,jfr3) + 360.0d0
      end if 
*-------------------------------------

      dega3(11,jfr3) = dega3(11,jfr3)*180.0d0/pi

*-------- correction of spline interpolation ---------
      if ((jfr3.ne.1).and.
     &           ((dega3(11,jfr3)-dega3(11,jfr3-1)).gt.180.0d0)) then
            dega3(11,jfr3) = dega3(11,jfr3) - 360.0d0
      else if ((jfr3.ne.1).and.
     &           ((dega3(11,jfr3)-dega3(11,jfr3-1)).lt.-180.0d0)) then
            dega3(11,jfr3) = dega3(11,jfr3) + 360.0d0
      end if 
*-------------------------------------

******Output to file****************************************
      open(unit=10,file='AnyBody/JointMotion/HipFlexionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega2(11,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega2(11,1)-dega2(11,jfr3)).gt.180.0d0)) then
             dega2(11,1) = dega2(11,1) - 360.0d0
            else if (((dega2(11,1)-dega2(11,jfr3)).lt.-180.0d0)) then
             dega2(11,1) = dega2(11,1) + 360.0d0
            end if 

            write(10,*) dega2(11,1)
           end if
*------------------------------------
        close(10)

      open(unit=10,file='AnyBody/JointMotion/HipAbductionLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega1(11,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega1(11,1)-dega1(11,jfr3)).gt.180.0d0)) then
             dega1(11,1) = dega1(11,1) - 360.0d0
            else if (((dega1(11,1)-dega1(11,jfr3)).lt.-180.0d0)) then
             dega1(11,1) = dega1(11,1) + 360.0d0
            end if

            write(10,*) dega1(11,1)
           end if
*------------------------------------
        close(10)

      open(unit=10
     &      ,file='AnyBody/JointMotion/HipExternalRotationLeft.any'
     &            ,status='unknown',access='append')
      write(10,*) dega3(11,jfr3),','
*-------- final frame data ----------
           if (jfr3.eq.frame) then

            if (((dega3(11,1)-dega3(11,jfr3)).gt.180.0d0)) then
             dega3(11,1) = dega3(11,1) - 360.0d0
            else if (((dega3(11,1)-dega3(11,jfr3)).lt.-180.0d0)) then
             dega3(11,1) = dega3(11,1) + 360.0d0
            end if

            write(10,*) dega3(11,1)
           end if
*------------------------------------
        close(10)


************** time for joint motion **********************

      open(unit=10,file='AnyBody/timeFrame.any'
     &            ,status='unknown',access='append')
      write(10,*) ttt*(jfr3-1.0d0)/frame,','
*-------- final frame data ----------
            if (jfr3.eq.frame) then
            write(10,*) ttt
           end if
*------------------------------------
        close(10)
*****************************************************


  515 continue

      return
      end
