calculs_skku_tan Subroutine

private subroutine calculs_skku_tan(bounds, lg, ssk, sku)

Note

Function to calculate the skewness and kurtosis of a tangent series

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:2) :: bounds

defines the function limits [-pi/2.(1-bounds(1)), +pi/2.(1-bounds(2)]
integer(kind=I4), intent(in) :: lg

vec size
real(kind=R8), intent(out) :: ssk

theoretical Ssk
real(kind=R8), intent(out) :: sku

theoretical Sku

Calls

proc~~calculs_skku_tan~~CallsGraph proc~calculs_skku_tan calculs_skku_tan proc~add_tang add_tang proc~calculs_skku_tan->proc~add_tang proc~tang tang proc~calculs_skku_tan->proc~tang proc~add_tang->proc~tang

Called by

proc~~calculs_skku_tan~~CalledByGraph proc~calculs_skku_tan calculs_skku_tan proc~calculs_skku_generique calculs_skku_generique proc~calculs_skku_generique->proc~calculs_skku_tan proc~fitness_skku_anal fitness_skku_anal proc~fitness_skku_anal->proc~calculs_skku_tan proc~cost_func_skku cost_func_skku proc~cost_func_skku->proc~fitness_skku_anal

Source Code

   subroutine calculs_skku_tan(bounds, lg, ssk, sku)
   !================================================================================================
   !<@note Function to calculate the skewness and kurtosis of a **tangent** series
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   real   (kind=R8), intent(in), dimension(1:2) :: bounds  !! *defines the function limits* [-pi/2.(1-bounds(1)), +pi/2.(1-bounds(2)]
   integer(kind=I4), intent(in)                 :: lg      !! *vec size*
   real   (kind=R8), intent(out)                :: ssk     !! *theoretical Ssk*
   real   (kind=R8), intent(out)                :: sku     !! *theoretical Sku*

      real(kind=R8)    :: xa, xb, mu, si, sk, ku, a, b
      real(kind=R8)    :: h, hh, b1, b2, alp, bet
      integer(kind=I4) :: i, ia, ib, deb, fin

      !---------------------------------------------
      ! WXMAXIMA file
      !---------------------------------------------
      !kill(all);
      !
      !f11(x):=tan(x)$
      !
      !assume(u<0.)$
      !assume(u>-%pi/2)$
      !assume(v>0.)$
      !assume(v<%pi/2)$
      !
      !I11:integrate(f11(x),x,u,v)$
      !I11:subst(-%pi/2+%pi/2*xa,u,I11)$
      !I11:subst(+%pi/2-%pi/2*xb,v,I11)$
      !I11:expand(trigsimp(I11));
      !
      !f21(x):=f11(x)-mu$
      !I21:integrate(expand(f21(x)^2),x,u,v)$
      !I21:subst(-%pi/2+%pi/2*xa,u,I21)$
      !I21:subst(+%pi/2-%pi/2*xb,v,I21)$
      !I21:expand(trigsimp(I21));
      !
      !f31(x):=f21(x)/si$
      !I31:integrate(f31(x)^3,x,u,v)$
      !I31:subst(-%pi/2+%pi/2*xa,u,I31)$
      !I31:subst(+%pi/2-%pi/2*xb,v,I31)$
      !I31:expand(trigsimp(I31));
      !
      !I41:integrate(f31(x)^4,x,u,v)$
      !I41:subst(-%pi/2+%pi/2*xa,u,I41)$
      !I41:subst(+%pi/2-%pi/2*xb,v,I41)$
      !I41:expand(trigsimp(I41));
      !---------------------------------------------

      ia = 256 ! ia and ib define the interval edges to be excluded ...
      ib = 256 ! ... because of high variations of the function.

      deb = 1  +ia ! start
      fin = lg -ib ! end

      a = bounds(1)
      b = bounds(2)

      hh = (2._R8-a-b)/(lg-1)
      h  = (PI_R8/2)*hh

      xa = a +ia*hh
      xb = b +ib*hh

      b1 = -PI_R8/2 *(UN-a)
      b2 = +PI_R8/2 *(UN-b)

      alp = -(b2-lg*b1)/(b2-b1)
      bet =      (lg-1)/(b2-b1)

      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      mu = log(sin((PI_R8*xa)/2.))-log(sin((PI_R8*xb)/2.))

      mu = (UN/h)*mu +add_tang(1, deb, fin, alp, bet, mu=0._R8, si=1._R8)
      do i = 1, ia-1
         mu = mu +tang(i*UN, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      do i = lg, lg -(ib -2), -1
         mu = mu +tang(i*UN, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      mu = mu/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      si = 2.*mu*log(sin((PI_R8*xb)/2.))+cos((PI_R8*xb)/2.)/sin((PI_R8*xb)/2.)-(PI_R8*mu**2.*xb)/2.+(PI_R8*xb)/2.                      &  !
          -2.*mu*log(sin((PI_R8*xa)/2.))+cos((PI_R8*xa)/2.)/sin((PI_R8*xa)/2.)-(PI_R8*mu**2.*xa)/2.+(PI_R8*xa)/2.+PI_R8*mu**2.-PI_R8      !

      si = (UN/h)*si +add_tang(2, deb, fin, alp, bet, mu, si=1._R8)
      do i = 1, ia-1
         si = si+ tang(i*UN, 2, alp, bet, mu, si=1._R8)
      enddo
      do i = lg, lg -(ib -2), -1
         si = si +tang(i*UN, 2, alp, bet, mu, si=1._R8)
      enddo
      si = si/lg
      si = sqrt(si)
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      sk =  log(sin((PI_R8*xb)/2.)**2.)/(2.*si**3.)                                 &  !
            -(3.*mu**2.*log(sin((PI_R8*xb)/2.)))/si**3.                             &  !
            -(3.*mu*cos((PI_R8*xb)/2.))/(si**3.*sin((PI_R8*xb)/2.))                 &  !
            +1/(2.*si**3.*sin((PI_R8*xb)/2.)**2.)+(PI_R8*mu**3.*xb)/(2.*si**3.)     &  !
            -(3.*PI_R8*mu*xb)/(2.*si**3.)                                           &  !
            -log(sin((PI_R8*xa)/2.)**2.)/(2.*si**3.)                                &  !
            +(3.*mu**2.*log(sin((PI_R8*xa)/2.)))/si**3.                             &  !
            -(3.*mu*cos((PI_R8*xa)/2.))/(si**3.*sin((PI_R8*xa)/2.))                 &  !
            -1/(2.*si**3.*sin((PI_R8*xa)/2.)**2.)+(PI_R8*mu**3.*xa)/(2.*si**3.)     &  !
            -(3.*PI_R8*mu*xa)/(2.*si**3.)-(PI_R8*mu**3.)/si**3.+(3.*PI_R8*mu)/si**3.   !

      sk = (UN/h)*sk +add_tang(3, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         sk = sk +tang(i*UN, 3, alp, bet, mu, si)
      enddo
      do i = lg, lg -(ib -2), -1
         sk = sk +tang(i*UN, 3, alp, bet, mu, si)
      enddo
      sk = sk/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ku = -(2.*mu*log(sin((PI_R8*xb)/2.)**2.))/si**4.+(4.*mu**3.*log(sin((PI_R8*xb)/2.)))/si**4.                                   &  !
           +(6.*mu**2.*cos((PI_R8*xb)/2.))/(si**4.*sin((PI_R8*xb)/2.))-(4.*cos((PI_R8*xb)/2.))/(3.*si**4.*sin((PI_R8*xb)/2.))       &  !
           -(2.*mu)/(si**4.*sin((PI_R8*xb)/2.)**2.)+cos((PI_R8*xb)/2.)/(3.*si**4.*sin((PI_R8*xb)/2.)**3.)                           &  !
           -(PI_R8*mu**4.*xb)/(2.*si**4.)+(3.*PI_R8*mu**2.*xb)/si**4.-(PI_R8*xb)/(2.*si**4.)                                        &  !
           +(2.*mu*log(sin((PI_R8*xa)/2.)**2.))/si**4.-(4.*mu**3.*log(sin((PI_R8*xa)/2.)))/si**4.                                   &  !
           +(6.*mu**2.*cos((PI_R8*xa)/2.))/(si**4.*sin((PI_R8*xa)/2.))-(4.*cos((PI_R8*xa)/2.))/(3.*si**4.*sin((PI_R8*xa)/2.))       &  !
           +(2.*mu)/(si**4.*sin((PI_R8*xa)/2.)**2.)+cos((PI_R8*xa)/2.)/(3.*si**4.*sin((PI_R8*xa)/2.)**3.)                           &  !
           -(PI_R8*mu**4.*xa)/(2.*si**4.)+(3.*PI_R8*mu**2.*xa)/si**4.-(PI_R8*xa)/(2.*si**4.)                                        &  !
           +(PI_R8*mu**4.)/si**4.-(6.*PI_R8*mu**2.)/si**4.+PI_R8/si**4.                                                                !

      ku = (UN/h)*ku +add_tang(4, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         ku = ku +tang(i*UN, 4, alp, bet, mu, si)
      enddo
      do i = lg, lg -(ib -2), -1
         ku = ku +tang(i*UN, 4, alp, bet, mu, si)
      enddo
      ku = ku/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ssk = sk
      sku = ku

   return
   endsubroutine calculs_skku_tan