DEEP2-redshift/contsubpoly.f at master · bjweiner/DEEP2-redshift · GitHub

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

c do continuum subtraction by subtracting a low order polynomial
c fit to the good data points.
c I added a sigma-clipping routine before the fit, which may not
c work great since the continuum varies, avsigclip or something
c would be better.

c do the subtraction in place, or not
c      subroutine contsubmed(spec,nsp,dwlog,wrad,specout)
      subroutine contsubpoly(wave,spec,serr,nsp,iord)

c      parameter (IORDER=15)

      real wave(nsp),spec(nsp),serr(nsp)
      real wscale(nsp),specclip(nsp)
      real wfit(nsp),ffit(nsp),efit(nsp)
      real a(iord),uu(nsp,iord),vv(iord,iord),ww(iord)
      real avals(iord)

      external polyfuncs

      include 'pcredshift.h'

      ncliprad=100
      clipsig=2.0
      call localsigclip(nsp,spec,specclip,ncliprad,clipsig)

      wrange = wave(nsp) - wave(1)
      do i=1,nsp
c rescale wfit to run between 0 and 1 for improved polynomial fitting
         wscale(i) = (wave(i) - wave(1)) / wrange
      end do

      npoly=iord
      nfit=0
      do i=1,nsp
         if(specclip(i) .gt. bad .and. serr(i) .gt. bad
     $        .and. serr(i) .lt. 1.0e6) then
            nfit=nfit+1
            wfit(nfit)=wscale(i)
            ffit(nfit)=specclip(i)
            efit(nfit)=serr(i)
         end if
      end do
      if (nfit .le. npoly) return

c do fit
      call bigsvdfit(wfit,ffit,efit,nfit,a,npoly,uu,vv,ww,nsp,iord,
     $     chisq,polyfuncs)

c evaluate and subtract fit
      do i=1,nsp
c         call polyfuncs(wave(i),avals,npoly)
         call polyfuncs(wscale(i),avals,npoly)
         fitval = 0.
         do j=1,npoly
            fitval = fitval + a(j)*avals(j)
         end do
         if (spec(i) .gt. bad) then
            spec(i) = spec(i) - fitval
         end if
      end do

      return
      end

cccccccccccccccccccccccccccccc
      subroutine polyfuncs(x,a,nord)
c IPOLY: 1 for simple polynomial: 1,x,x^2,...
      parameter(IPOLY=1)

      real a(nord)

      if (IPOLY .eq. 1) then
         a(1) = 1.0
         do i=2,nord
            a(i) = a(i-1)*x
         end do
      end if

      return
      end