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?

一日の中央値を求めるプログラム

使い方

$ gfortran -o medi1day sac_median_1day.f95
$ ./medi1day ファイル名

sac_median_1day.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

  if(iargc()>0) then
    call getarg (1, file1) !コマンドラインで引数
 endif

 max=10000000
 allocate (x(max),y(max))
 call arsac(file1,b,a,khd,x,max)
  npts = a(10)
  dt = b(1)
  call median(x, med, 8640000)
  med = med*1e-9 !nm/sからm/sへ変換
 write(*,*)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,8640001)

  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?