Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
0
Help us understand the problem. What are the problem?

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

使い方

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
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
0
Help us understand the problem. What are the problem?