/
Correl.f
executable file
·140 lines (111 loc) · 3.91 KB
/
Correl.f
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
subroutine correl
c******************************************************************************
c This routine cross-correlates arrays x and y using any part of the
c arrays, as long as there are enough data points in each array.
c The result is printed out in terms of the shift in the abscissa
c of the y array needed to align it with the x array.
c******************************************************************************
include 'Atmos.com'
include 'Linex.com'
include 'Pstuff.com'
include 'Equivs.com'
include 'Plotval.com'
real*4 xgood(100000), ygood(100000), zgood(100000)
real*8 spx(1000), spy(1000)
real*4 q
real*8 xinterp(21), yinterp(21)
wavemin = amax1(xsyn(1),xobs(1))
wavemax = amin1(xsyn(kount),xobs(lount))
wavecenter = (wavemax + wavemin)/2.
c*****find the array subscripts for the minimum and maximum
c wavelengths in the synthetic spectrum.
do i=1,kount
if (xsyn(i) .ge. wavemin) then
ilosyn = i
go to 5
endif
enddo
5 do i=ilosyn,kount
if (xsyn(i) .gt. wavemax) then
ihisyn = i - 1
go to 10
endif
enddo
ihisyn = kount
c*****dump the synthetic spectrum wavelengths and fluxes into working arrays
10 itotsyn = ihisyn - ilosyn + 1
do i=1, itotsyn
xgood(i) = xsyn(ilosyn-1+i)
ygood(i) = chunk(ilosyn-1+i,1)
enddo
c*****find the array subscript for the minimum wavelength in the
c observed spectrum
do j=1,lount
if (xobs(j) .ge. xgood(1)) then
jloobs = j
go to 15
endif
enddo
c*****using a 3-point Lagrangian, make an interpolated observed
c spectrum that matches the wavelength step size of the
c synthetic spectrum
15 j = jloobs
do i=1,itotsyn
if (xgood(i) .gt. (xobs(j+1)+xobs(j))/2.) j = j + 1
q = xgood(i)
zgood(i) = yobs(j-1)*(q-xobs(j))*(q-xobs(j+1))/
. ((xobs(j-1)-xobs(j))*(xobs(j-1)-xobs(j+1))) +
. yobs(j)*(q-xobs(j-1))*(q-xobs(j+1))/
. ((xobs(j)-xobs(j-1))*(xobs(j)-xobs(j+1))) +
. yobs(j+1)*(q-xobs(j-1))*(q-xobs(j))/
. ((xobs(j+1)-xobs(j-1))*(xobs(j+1)-xobs(j)))
enddo
c*****next, EITHER do a single cross-correlation
call crosscorr (ishift,spy(1),ygood,zgood,itotsyn)
c*****next, do the cross-correlation:
imax = 2*maxshift + 1
do i=1,imax
ishift = i - 1 - maxshift
spx(i) = dble(ishift)
call crosscorr (ishift,spy(i),ygood,zgood,itotsyn)
enddo
c*****interpolate to find the maximum of the correlation function
corrmax = spy(1)
mc = 1
do i=2,imax
if (spy(i) .gt. corrmax) then
corrmax = spy(i)
mc = i
endif
enddo
if (mc.eq.1 .or. mc.eq.imax) then
deltawave = 0.
return
endif
xinterp(1) = spx(mc-1)
yinterp(1) = spy(mc-1)
do j=2,21
q = spx(mc-1) + (j-1)/10.
xinterp(j) = q
yinterp(j) = spy(mc-1)*(q-spx(mc))*(q-spx(mc+1))/
. ((spx(mc-1)-spx(mc))*(spx(mc-1)-spx(mc+1))) +
. spy(mc)*(q-spx(mc-1))*(q-spx(mc+1))/
. ((spx(mc)-spx(mc-1))*(spx(mc)-spx(mc+1))) +
. spy(mc+1)*(q-spx(mc-1))*(q-spx(mc))/
. ((spx(mc+1)-spx(mc-1))*(spx(mc+1)-spx(mc)))
enddo
jmax = 1
corrmax = yinterp(1)
do j=2,21
if (yinterp(j) .gt. corrmax) then
jmax = j
corrmax = yinterp(j)
shiftmax = xinterp(j)
endif
enddo
c*****translate this into a correction to the input velocity shift
deltawave = -step*shiftmax
deltavel = -3.0d5*step*shiftmax/wavecenter
veladd = veladd + deltavel
return
end