LoginSignup
0
0

More than 1 year has passed since last update.

3年分の振幅計算(1時間の中央値)

Posted at

使い方

N.MWDH観測点の2-5Hz,15時〜一時間の中央値の場合
./allmedian.sh _N.MWDH_U_bp2-5_sqr 15

allmedian.sh
#!/bin/bash

#./allmedian.sh _N.MMOH_U_bp2-5_sqr 15

for yr in {19..20}
do

 for month in {1..9}
 do
  for day in {1..9}
  do
  ./medi2 ${yr}0${month}/${yr}0${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done

  for day in {10..31}
  do
  ./medi2 ${yr}0${month}/${yr}0${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done
 done

 for month in {10..12} 
 do
  for day in {1..9}
  do
  ./medi2 ${yr}${month}/${yr}${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done

  for day in {10..31}
  do
  ./medi2 ${yr}${month}/${yr}${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done
 done

done
########
yr=21

 for month in {1..5}
 do
  for day in {1..9}
  do
  ./medi2 ${yr}0${month}/${yr}0${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done

  for day in {10..31}
  do
  ./medi2 ${yr}0${month}/${yr}0${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
  done
 done

# for month in {10..12} 
#
#  for day in {1..9}
#  ./medi2 ${yr}${month}/${yr}${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
#  done
#
#  for day in {10..31}
#  ./medi2 ${yr}${month}/${yr}${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
#  done
# done

スクリプト中に含まれるmedi2は以下のプログラムから

sac_median_2.f95
program main
 implicit none
 integer i, a(40), max, npts, hour
 real b(70),c(2), dt, med
 real, allocatable :: x(:), y(:)
 character :: file1*80, arg*10, khd(48)*4, name*10

  if(iargc()>0) then
    call getarg (1, file1) !コマンドラインで引数!
    call getarg (2, arg) !コマンドラインで引数所得
    read(arg,*) hour
 endif

 max=10000000
 allocate (x(max),y(max))
 call arsac(file1,b,a,khd,x,max)
  npts = a(10)
  dt = b(1)
  do i = hour*360000+1,(hour+1)*360000+1
     y(i-hour*360000) = x(i)
  enddo
  call median(y, med, 360000)
  med = med*1e-9 !nm/sからm/sへ変換
!  write(*,*)med
  name = '20'//file1(6:11)
  write(*,*)name, med
end program main
!########################################################
subroutine arsac(fn,fhd,nhd,khd,y,maxy)
  real fhd(70),y(maxy)
  integer nhd(40),i,i1,i2,i3,npts
  character khd(48)*4,fn*80

  open(1,file=fn,access='direct',recl=632,status='old')
  read(1,rec=1)(fhd(i1),i1=1,70),(nhd(i2),i2=1,40),(khd(i3),i3=1,48)
  npts=nhd(10)
  close(1)
  open(1,file=fn,access='direct',recl=4,status='old')
  do i=1,npts
!     read(1,rec=i+158,end=100)y(i)     Mar.1 H.Aoyama
     read(1,rec=i+158)y(i)
  enddo
! 100 continue
  close(1)
  return
end subroutine arsac
!#####################################################
subroutine median(a, med, n)
  implicit none
  integer,intent(in):: n
  real,intent(inout):: a(n)
  real,intent(inout):: med
  integer i, m
  real am

!バブルソート(遅い)
!  do i = 1, n-1 !sort a
!     am = minval(a(i+1:n))
!     m  = minloc(a(i+1:n), 1) + i
!     if (a(i) > am) then
!        a(m) = a(i)
!        a(i) = am
!     endif
!  enddo

  call quicksort(a,1,360001)

  if(mod(n, 2) /= 0) then !median of a(:)
     med = a((n-1)/2+1)
  else
     med = 0.5*(a(n/2)+a(n/2+1))
  endif
end subroutine median

!##############################################
recursive subroutine quicksort(a, first, last)
  implicit none
  real  a(*)
  integer first, last
  integer i, j
  real  x, t

  x = a( (first+last) / 2 )
  i = first
  j = last
  do
     do while (a(i) < x)
        i=i+1
     end do
     do while (x < a(j))
        j=j-1
     end do
     if (i >= j) exit
     t = a(i);  a(i) = a(j);  a(j) = t
     i=i+1
     j=j-1
  end do
  if (first < i - 1) call quicksort(a, first, i - 1)
  if (j + 1 < last)  call quicksort(a, j + 1, last)
end subroutine quicksort
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0