!----------------- 平滑子程序 --------------------------- !
!本文件共有两个子程序:                                    !
!       1.sub_ph1(n,ki,x,y):线性平滑;                     !
!       2.sub_ph2(n,ki,x,y):二项式平滑;                   !
!要求输入的有:                                            !
!       1.x(n):长度为n的序列;                             !
!       2.ki:平滑步长;                                    !
!输出的有:                                                !
!       1.y(n):平滑后的序列;                              !
!------------------------------------- 程正泉 2000.5 -----!

!=========================================================!
subroutine sub_ph1(n,ki,x,y)                              !
implicit none                                             !
integer,intent(in)::n,ki                                  !
real,dimension(n),intent(in)::x                           !
real,dimension(n),intent(out)::y                          !
integer::i,j                                              !
real::tmp                                                 !
do i=1,n                                                  !
   if((i >= int(ki/2)) .and. (i <= n-int(ki/2)))then      !
     tmp=0.0                                              !
     do j=i-int(ki/2),i+int(ki/2)                         !
        tmp=tmp+x(j)                                      !
     enddo                                                !
     y(i)=tmp/real(ki)                                    !
   else                                                   !
     y(i)=0.0                                             !
   endif                                                  !
enddo                                                     !
end                                                       !
!=========================================================!

!=========================================================!
subroutine sub_ph2(n,ki,x,y)                              !
implicit none                                             !
integer,intent(in)::n,ki                                  !
real,dimension(n),intent(in)::x                           !
real,dimension(n),intent(out)::y                          !
real,allocatable,dimension(:)::h                          !
integer::i,j                                              !
real::temp                                                !
                                                          !
allocate(h(-ki/2:ki/2))                                   !
call coef(h,ki)                                           !
do i=1,n                                                  !
  if(i<=ki/2)then                                         !
    y(i)=0.0; temp=0.0                                    !
    do j=1-i,ki/2                                         !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  elseif((i > ki/2).and.(i < n-Ki/2+1))then               !
    y(i)=0.0; temp=0.0                                    !
    do j=-ki/2,ki/2                                       !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  else                                                    !
    y(i)=0.0; temp=0.0                                    !
    do j=-ki/2,n-i                                        !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  endif                                                   !
enddo                                                     !
deallocate(h)                                             !
end                                                       !
!     ---     ---     ---     ---    ---    ---    ---    !
subroutine coef(h,k)                                      !
!二项式的系数,参考丁裕国教授所编著的气象...               !
implicit none                                             !
integer,intent(in)::k                                     !
real,dimension(-k/2:k/2),intent(out)::h                   !
integer,external::ppp                                     !
integer::i                                                !
do i=-k/2,k/2                                             !
  h(i)=ppp(k-1)*1.0/(ppp(i+k/2)*ppp(k-1-i-k/2))           !
enddo                                                     !
end                                                       !
!     ---     ---     ---    ---     ---     ---    ---   !
!求n的阶乘(n!)                                            !
integer function ppp(n)                                   !
implicit none                                             !
integer,intent(in)::n                                     !
integer::i                                                !
if(n == 0)then                                            !
  ppp=1                                                   !
else                                                      !
  ppp=1                                                   !
  do i=1,n                                                !
    ppp=ppp*i                                             !
  enddo                                                   !
endif                                                     !
end                                                       !
!=========================================================!