1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
!Copyright (C) 2013- Roger Wesson
!Free under the terms of the GNU General Public License v3
program cont
use mod_quicksort
implicit none
character(len=50) :: filename, windowch
character(len=1) :: z
integer :: i, io, spectrumlength, begin, end, window
type spectrum
real(kind=dp) :: wavelength
real(kind=dp) :: flux
real(kind=dp) :: difference
real(kind=dp) :: continuum
end type
type(spectrum), dimension(:), allocatable :: input
type(spectrum), dimension(:), allocatable :: chunk
real(kind=dp), dimension(:), allocatable :: temporary
real(kind=dp) :: meanflux, meandifference
window=25
!continuum points have the lowest flux and represent low frequency information so will have low gradients to adjacent points.
!so, first select lowest half of the fluxes
!then select lowest half of differences to points either side.
call get_command_argument(1,filename)
if (filename .eq. "") then
print *,"file?"
read(5,*) filename
endif
if (iargc() .eq. 2) then
call get_command_argument(2,windowch)
read (windowch,"(I3)") window
endif
I = 0
OPEN(199, file=filename, iostat=IO, status='old')
DO WHILE (IO >= 0)
READ(199,*,end=112) z
I = I + 1
END DO
112 spectrumlength=I
!then allocate and read
allocate(input(spectrumlength))
allocate(chunk(2*window+1))
allocate(temporary(2*window+1))
REWIND (199)
DO I=1,spectrumlength
READ(199,*) input(i)%wavelength,input(i)%flux
input(i)%continuum = 0.d0
END DO
CLOSE(199)
!now get first differences
do i=2,spectrumlength-1
input(i)%difference = abs(input(i)%flux - input(i-1)%flux) + abs(input(i+1)%flux - input(i)%flux)
end do
!in 100 unit windows, take continuum as average of points with less than mean 1st difference and less than mean flux.
do i=1,spectrumlength
if (i-window .lt. 1) then
begin=1
else
begin=i-window
endif
if (i+window .gt. spectrumlength) then
end=spectrumlength
else
end=i+window
endif
chunk = input(begin:end)
meandifference=sum(chunk%difference)/size(chunk)
meanflux=sum(chunk%flux)/size(chunk)
where (chunk%difference .gt. (meandifference/2.))
chunk%continuum=0.d0
elsewhere
chunk%continuum=chunk%flux
endwhere
input(i)%continuum=sum(chunk%continuum)/count(chunk%continuum.gt.0.d0)
enddo
open(123, file="continuum.dat")
do i=1,spectrumlength
write (123,*) input(i)%wavelength,input(i)%flux,input(i)%continuum
enddo
end program cont
|
