c
c##########################################################################
c Transit search with robust TFA+FOUR filtering for stellar variability
c and outlier/transit preservation
c##########################################################################
c
c DATE: - September 07, 2020
c
c INPUT: - tran_k2_v0.par :: Input parameter file
c - filename.lis :: List of stars to be analyzed with complete paths
c to their time series
c - ref.pos :: Reference position/magnitude data file of all
c stars available
c - template.lis :: List of templates & paths to their time series
c
c OUTPUT: - snr.dat :: Result of the analysis on the stars analyzed
c [this is an "appended" file]
c - spec/*.sp :: Frequency spectrum of the star analyzed
c - flc/*.lc :: Filtered [non-reconstructed] LC of the star analyzed
c
c - sp.dat :: Frequency spectrum of the star analyzed
c - four_ord.dat :: Diagnostic parameters for the optimum Fourier order
c - test_par.dat :: Test signal parameters ( produced only if itest > 0 )
c [this is an "appended" file]
c - ts_4.dat :: Time series output [after TFA+FOUR filtering]
c - outlier.dat :: Outlying points as given by the routine "outlier1"
c - pre_bls.dat :: Data AFTER AR edge correction for Gibbs phenomenon
c (this is the time series used by the BLS search)
c - sp##.dat :: BLS spectr. after the ##-th prewhitening [no reconstruction]
c - lc##.dat :: LC after the ##-th prewhitening [no reconstruction]
c - pa##.dat :: Folding params for the ##-th prewhitening [no reconstruction]
c - rec_tfa.dat :: LC after the SINGLE transit reconstruction step
c - multi_bls.dat :: BLS parameters after the FULL prewhitening cycle
c - multi_par.dat :: Parameters of the significant components derived by the
c iterative reconstruction
c - tr#.lc :: Individual/reconstructed transit components and the
c corresponding fits with folding phase values
c - tr_all.lc :: Reconstructed transit LC and the fit with all components
c - tr#_fit.lc :: Densely sampled transit fit to the #-th component
c
c
c PLOTS: - cand_4_plot.png:: 4-panel diagnostic plot for the main transit component of
c a given candidate
c - mult_plot.png :: BLS spectra and zoomed/folded LCs for each transit
c component of a given candidate
c - plots/*.png :: The above files with star IDs
c
c
c NOTES: - This is the "cleaned-up" version of the code "tran_k2" that was used
c to produce all results presented in the paper on this code and the
c survey of K2/C05 (submitted on 2020.06.23).
c
c - Not all output files are generated after each run. For example, if
c isurvey=1, then the only output file produced is *snr.dat*
c [and, if itest > 0, then the test signal parameters]
c
c - The same type of data product might appear in different files at
c various stages of filtering. This may be useful in checking the
c consistency of these data products.
c
c - Plots are ALWAYS generated, except if isurvey=0. They are saved under
c target-specific names only if ispl=1
c
c - Plots are generated by using GNUPLOT. Files required:
c * fold_k2_rr.f
c * cand_4_plot.sh
c * mult_plot.sh
c * mult_#.sh - with #=1,2,3,4,5
c
c - Maximum number of time series data points = 30333
c Maximum number of frequency spectrum values = 155555
c
c - There is a seemingly harmless error message at the end of each run
c that we could not cure:
c "The following floating-point exceptions are signalling:
c IEEE_DIVIDE_BY_ZERO"
c
c=================================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fname
character*9 sname,ref_id(30333)
character*3 camp
character*70 command
c
dimension t(30333),x(30333),we(30333),x_rec(30333)
dimension xf(30333),y(30333),xw(30333),u(30333),v(30333)
dimension told(30333),xold(30333),xsyn(30333),x_int(30333)
dimension ref_mag(30333),ref_ra(30333),ref_de(30333)
dimension xft(30333),fourts(30333),tfats(30333)
dimension xf1(30333),xf2(30333),xp1(30333),xp2(30333)
dimension tp1(30333),tp2(30333),xfa(30333),z(30333),xtr(30333)
dimension trsyn(30333),x_1(30333),tfa_r(30333),four_r(30333)
dimension ft(6,30333),xt(6,30333)
dimension fw(100),dw(100),snrw(100),fr(100)
dimension td(100),qt(100),ri(100),tc(100)
dimension p(155555)
dimension intr(30333)
c
common /tt/ ttemp(30333),ntemp
common /rc/ m3_rc
c
c=================================================================================
c
c... WIRED-IN PARAMETERS
c
tc0 = 57000.0d0
ntw = 6
pfac = 1.05d0
ext = 0.05d0
mar = 900
iar = 1
npol = 5
iew = 0
iwa = 0
ndsp = 6
nb_sp = 2000
iso = 1
isel = 1
asig2 = 3.0d0
snr_min = 6.0d0
ising = 0
dsp_min = 10.0d0
ksm = 1
c
c tc0 :: Transit epoch for the test signal
c ntw :: Number of prewhitening cycles for the multiple transit search
c pfac :: The fundametal period of the Fourier component is
c pfac*(t_max-t_min), to avoid (or decrease) the Gibbs
c phenomenon at the edges.
c ext :: AR prediction is performed on both sides for N*ext
c data points at the edges
c mar :: Order of the AR process
c iar :: 0 - do not use AR process to decrease Gibbs oscillations
c 1 - use AR process to decrease Gibbs oscillations
c npol :: Order of the polynomial for filtering out trends from the
c AR predictions
c iew :: 1 - use equal weights in "one_rec"
c 0 - use weights in "one_rec" as given in that routine, but
c do not employ "exact zero weight" windows.
c 2 - as iew=0, but DO EMPLOY "exact zero weight" windows.
c iwa :: 0 - do not employ exact ZERO weights
c 1 - employ exact ZERO weights [when appropriate - see "slide_rms"]
c NOTE: This parameter is used ONLY in the basic signal processing
c ndsp :: Order of the polynomial fit to the BLS spectrum
c nb_sp :: Number of frequency bins for the economically sampled frequency
c spectra
c iso :: 0/1 - no/yes single outlier correction
c isel :: 0 - Take always the type of edge correction given by "iar"
c iar=0 -- we have FOUR edge fit
c iar=1 -- we have AR edge fit
c 1 - Compare the FOUR and AR approximations and take the one
c with the smallest residual robust scatter.
c iar=0 -- we have FOUR edge fit
c iar=1 -- we have either FOUR or AR edge fit, depending
c on the size of the residual robust scatter
c asig2 :: Sigma clipping parameter. Employed ONLY on the reconstructed,
c multi-component time series (once all the components are found.
c snr_min:: Minimum spectrum SNR to declare a transit component detected
c in the multi transit reconstruction routine "multi_rec".
c ising :: 1 - Perform single transit test
c 0 - Skip single transit test
c dsp_min:: Minimum DSP to declare a single transit event as such.
c ksm :: 1- Employ reconstruction in the multi-tranist parameter
c determination
c 0- Do NOT employ reconstruction in the multi-tranist parameter
c determination [this option DOES not work yet]
c
c=================================================================================
c
c... READ IN BASIC PARAMETERS
c
call parin(isurvey,jsm,itemp,nstar,fmin1,fmax1,nf1,
&iwhite,nb,qmi,qma,itest,tfreq0,tfreq1,dep0,dep1,t14_0,
&t14_1,pht0,pht1,frin0,frin1,amootv,phootv,frootv,iflare,
&igen,std,fdep,iseed,issp,isflc,ispl,iplot,asig,ic_flare,
&ibg,ixy,m_four,m_opt)
c
c.................................... Set certain parameters
iseed0= iseed
iseed = 2000*iseed+1
totdf = fmax1-fmin1
fmin11= fmin1
nbb = 2000
qtr0 = t14_0*tfreq0
qtr1 = t14_1*tfreq1
dep00 = dep0
dep10 = dep1
rfac = 0
c
c.................................... Change some parameters
if(iplot.eq.0) ispl=0
if(isurvey.eq.0) go to 11
issp=0
isflc=0
iwhite=0
11 continue
c
c... Read in *ref.pos*
c
call read_ref(n_ref,ref_id,ref_ra,ref_de,ref_mag)
c
c... Create directories "spec", "flc" and "plots" if:
c a) they do not exisist
c b) Any of these parameters is greater than zero:
c isflc, issp, ispl
c
command='if [ ! -d spec ]; then mkdir spec; fi'
if(issp.gt.0) call SYSTEM(command)
command='if [ ! -d flc ]; then mkdir flc; fi'
if(isflc.gt.0) call SYSTEM(command)
command='if [ ! -d plots ]; then mkdir plots; fi'
if(ispl.gt.0) call SYSTEM(command)
c
open(14,file='snr.dat',access='append')
if(itest.gt.0) open(12,file='test_par.dat',access='append')
write(14,2)
2 format('# SNAME N TOT AVMAG SIG_FINAL',
&' F0 DEPTH SNR SPD DSP QTRAN TCEN',
&' NEV NDIT SIG_FOUR SIG_TFA SNRT FRET',
&' RING ITOC F03 SNR03 F35 SNR35 M_FOUR')
c
c
c**************************************
c START CYCLE OF ANALYSIS *
c**************************************
c
do 1 nseq=1,nstar
c
write(*,6) nseq,nstar
6 format('*************************** CYCLE ',2(1x,i6))
c
call cpu_time(t_start)
c
if(itest.eq.0) go to 15
if(nstar.eq.1) go to 15
rfac=2.0d0*arand(iseed)-1.0d0
dep0=dep00*(1.0d0+rfac*fdep)
dep1=dep10*(1.0d0+rfac*fdep)
15 continue
c
c... Read in the first file item of 'filename.lis'
call file_in(fname,sname,nfile,camp)
c
c... Match object with its average magnitude
call match_obj(sname,n_ref,ref_id,ref_mag,avmag)
c
c... Read in data (time series, etc.) and perform TFA+FOUR filtering
call datain_tfa_four(isurvey,nseq,itest,iseed,m_four,
&itemp,fname,sname,n,t,x,nvold,told,xold,
&intemp,avmag,frin0,frin1,pht1,tc0,frootv,amootv,
&phootv,qtr0,qtr1,pht0,igen,std,xsyn,trsyn,
&x_int,camp,we,iflare,ibg,ixy,sig,pfac,ext,
&iwa,fourts,tfats,tfreq0,dep0,tfreq1,dep1,m_opt)
c
c... Brief summary:
c - {n, t(i),x_int(i)} :: ORIGINAL (interpolated to {ttemp}) time series
c - {tfats} :: TFA vector
c - {fourts} :: FOURIER vector
c - {x} :: Filtered time series: {x_int}-{tfats}-{fourts}
c
if(intemp.eq.1) write(*,151) sname
151 format(/,'>>>>>>>>> Target ',a10,
&' was a member of the template set!',/)
c
do 9 i=1,n
xfa(i)=fourts(i)
xf(i)=fourts(i)
9 continue
npr=0
if(iar.eq.0) go to 150
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c AUTOREGRESSIVE EDGE CORRECTION FOR THE GIBBS OSCILLATIONS
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c... Generate noisy truncated signal from the Fourier fit
c ... We use UNIFORM noise distribution [to avoid large
c perturbations in the synthetic Fourier signal
c [ this is why we have that sqrt(3) factor ].
c
npr=idint(dfloat(n)*ext)
sigf=dsqrt(3.0d0)*sig
do 8 i=1,n
z(i)=x_int(i)-tfats(i)
y(i)=xf(i)+sigf*(2.0d0*arand(iseed)-1.0d0)
8 continue
c
c.................................................... Cut edges for AR modeling
call cut_ts(npr,n,t,y,n_ar,u,v)
c
c.................................................... AR prediction at the edges
call LR_pred(mar,npr,n_ar,u,v,xf1,xf2,xp1,xp2,
&tp1,tp2,xw,sig)
c
c.................................................... Full - AR+FOUR+AR - time series {xfa}
call ar_four_ar(npr,n,t,xf,tp1,xp1,tp2,xp2,xfa)
c
c.................................................... Robust fit of low-order polynomials
c to the AR extrapolations at the edges
call end_pol(npol,npr,n,t,z,xfa)
c
c.................................................... Compare FOUR and AR and select the better one
if(isel.eq.1) call ar_vs_four(npr,n,z,xf,xfa,s_xf1,s_xfa1,
&s_xf2,s_xfa2,iflag1,iflag2)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
150 continue
c
c The time series below, {t(i),x(i)}, is corrected for
c systematics {tfats} and stellar variability {xfa}.
c
c HOWEVER: it is still *incomplete/biased*, if there are
c other signals, i.e., transits. Correction for these
c additional components will be done by the
c reconstruction routines, e.g., by *one_rec*.
c-----------------------------------------------------------------
do 10 i=1,n
x(i) = x_int(i)-xfa(i)-tfats(i)
10 continue
c
c........................................ Omitting single outliers
if(iso.eq.1) call outlier_1g(isurvey,1,n,t,x,nout)
c
c........................................ Flare correction
if(ic_flare.eq.1) call kill_flare(1,n,x,nflare)
c
c........................................ Sigma clipping, zero-average, save PRE-BLS
call sigma_clipping(1,asig,n,x,nclip)
c
c........................................ Save PRE-BLS data and set {x_rec}={x}
if(isurvey.ne.1) call save_pre_bls(isflc,sname,n,t,x,xfa,
&tfats,we,x_rec)
c
c........................................ Check efficacy of systematics and
c stellar variability removal
call check_sys_four(n,t,x,f03,snr03,f35,snr35)
c
c........................................ Compute RMS of the FOUR and TFA
c corrections
call four_vs_tfa(n,tfats,xfa,sig_four,sig_tfa)
c
c........................................ Save the AR correction
if(isurvey.eq.0) call ts_ar_dat(npr,n,t,x_int,xf,xfa,x)
c
c................................. Compute total time span and frequency range
call minmax(n,t,tmin,tmax)
tot = tmax-tmin
fmin0 = 1.001d0/tot
fmin1 = fmin11
if(fmin1.lt.fmin0) fmin1=fmin0
fmax1 = fmin1+totdf
df1 = (fmax1-fmin1)/dfloat(nf1-1)
c
c................................. Search for single transit event
nsing=0
if(ising.eq.1) call single_tran(n,t,x,tmin,tmax,dsp_min,
&tce_1,dep_1,qtr_1,qin_1,x_1,nev_1,ndit_1,itoc_1,dsp_1,nsing)
c
if(nsing.eq.0) go to 50
write(*,52) dep_1,qtr_1,qin_1,tce_1,nev_1,ndit_1,itoc_1,dsp_1
52 format('Single transit parameters:',/,
&'Depth = ',f13.6,/,
&'Q_tran = ',f13.6,/,
&'Q_ingr = ',f13.6,/,
&'T_cen = ',f13.5,/,
&'N_event= ',i13,/,
&'N_intr = ',i13,/,
&'ITOC = ',i13,/,
&'DSP = ',f13.1)
50 continue
c
c................................. Standard BLS search
call eebls(n,t,x,u,v,nf1,fmin1,df1,nb,qmi,qma,
&p,bper,bpow,depth,qtran,in1,in2)
c
c................................. Filter out low-frequency trend from the frequency spectrum
call sppolfit_w(ndsp,nf1,p)
c
c................................. Save frequency spectrum
call save_sp(nb_sp,nf1,fmin1,df1,p)
if(issp.eq.1) call save_spec(sname,nb_sp,nf1,fmin1,df1,p)
c
c................................. Compute frequency at the highest peak and SNR
call spsnr1(nf1,p,fmin1,df1,frmax0,pmax0,sde0,spd0)
c
if(nsing.eq.0) go to 54
p0=tmax-tmin
f0=1.0d0/p0
frmax0=f0
go to 51
54 continue
c
c................................. Compute ACCURATE period "p0", check 1st (sub)harmonic
p10=0.5d0/frmax0
pac=p10
dev=1.0d30
do 30 j=1,3
p1=p10*(2**(j-1))
call ac_per(p1,df1,n,t,x,nbb,qmi,qma,pac1,sig)
if(sig.gt.dev) go to 30
dev=sig
pac=pac1
30 continue
p0=pac
f0=1.0d0/p0
c
51 continue
c
c................................. Light curve reconstruction and transit parameter determination
nout=0
nflare=0
nclip=0
call one_rec(isurvey,iew,jsm,p0,n,t,x_int,x,nbb,qmi,qma,
&m_four,iso,intr,nintr,depth,ring,qtran,tcen,xft,x_rec,we,
&sig,nout,tfa_r,four_r)
c
c................................. Flare correction
if(ic_flare.eq.1) call kill_flare_r(2,n,x_rec,xft,nflare)
c
c... NOTE: We do not have relative sigma clipping in respect of
c {xft} - sigma_clipping_r(2,asig2,n,x_rec,xft,nclip) -,
c because this would clip also possible other components.
c
c................................. Compute DSP, SNRT, FRET
call dsp_ok(n,t,x_rec,p0,tcen,qtran,dsp,snrt,fret)
c
c................................. Compute unbiased estimate of the standard deviation of
c the residuals around the best-fitting TRANSIT model
do 550 i=1,n
y(i)=x_rec(i)-xft(i)
550 continue
call astat(y,n,aver,sigma)
n_red=n-4-m3_rc-nout-nflare-nclip
sig_final=sigma*dsqrt(dfloat(n)/dfloat(n_red))
c
c................................. Save filtered/reconstructed time series
call save_rec(n,t,x_rec,xft,xsyn,trsyn,
&tfa_r,four_r,qtran)
c
if(isurvey.eq.1) go to 14
c
c................................. Save peak frequency parameters in "fold.par"
open(24,file='fold.par')
nper=1
write(24,40) sname,nper,f0,tcen
40 format(a9,/,i2,/,f11.7,/,f11.4)
close(24)
c
c................................. Create and save basic diagnostic plot
if(iplot.ne.0) call SYSTEM('gnuplot cand_4_plot.sh')
command=trim('cp cand_4_plot.png plots/cand_4_'//sname//'.png')
if(ispl.eq.1) call SYSTEM(command)
c
14 continue
c
c................................. Compute number of transits and the number of the observed
c data points in them
call events(n,t,f0,tcen,qtran,nev,ndit,itoc)
c
c................................. Print out result of the analysis of the TFA+FOUR+AR-filtered data
write(*,3) sname,n,tot,avmag,sig_final,f0,depth,sde0,spd0,
&dsp,qtran
write(14,3) sname,n,tot,avmag,sig_final,f0,depth,sde0,spd0,
&dsp,qtran,tcen,nev,ndit,sig_four,sig_tfa,snrt,fret,ring,itoc,
&f03,snr03,f35,snr35,m_four
3 format(a9,2x,i5,1x,f5.0,1x,f8.5,1x,f9.6,1x,f12.8,1x,f9.6,
&1x,f4.1,1x,f5.3,1x,f4.1,1x,f6.4,1x,f14.7,2(1x,i5),2(1x,f8.6),
&1x,f6.2,1x,f7.2,1x,f7.4,1x,i5,2(1x,f10.6,1x,f5.1),5x,i3)
c
c................................. Print out test data
if(itest.gt.0) write(12,20) sname,tfreq0,dep0,t14_0,pht0,frin0,
&tfreq1,dep1,t14_1,pht1,frin1,amootv,phootv,frootv,iflare,igen,
&std,1.0d0+rfac*fdep,iseed0
20 format(a9,2(2x,f12.8,1x,f9.6,1x,f5.3,1x,f5.3,1x,f5.3),
&1x,f9.6,1x,f5.3,1x,f12.8,2(1x,i2),1x,f9.6,1x,f5.3,1x,i6)
c
if(iwhite.eq.0) go to 5
c
c#################################################################
c SUCCESSIVE PREWHITENING BY BLS SIGNALS
c
call bls_white(n,t,x,nf1,fmin1,df1,nb,qmi,qma,
&nsing,p0,sname,nb_sp,ndsp,nbb,ntw,fw,dw,snrw)
c
c................................. Reconstruction of the significant transit components
call multi_rec(ksm,n,t,x_int,x,nbb,
&qmi,qma,ntw,fw,dw,snrw,snr_min,asig2,
&ntc,xtr,xft,ft,xt,we,fr,td,qt,ri,tc,sig,nclip)
c
c................................. Save reconstructed LCs and transit parameters
call save_m_rec(sname,ntc,n,t,xtr,xft,ft,xt,
&fr,td,qt,ri,tc)
c
c#################################################################
c
c................................. Create and save plots for the significant transit components
if(iplot.ne.0) call SYSTEM('./mult_plot.sh')
command=trim('cp mult_plot.png plots/mult_'//sname//'.png')
if(ispl.eq.1) call SYSTEM(command)
c
5 continue
c
c................................. Print out execution time
call cpu_time(t_finish)
write(*,4) t_finish-t_start
4 format(/,'EXECUTION TIME [s] = ',f7.2,/)
c
1 continue
c
close(14)
if(itest.gt.0) close(12)
c
stop
end
c
c
subroutine sigma_clipping(ipr,asig,n,x,nclip)
c--------------------------------------------------------------------------
c - Performs iterative sigma-clipping of {n,x} at asig*sigma level
c - Upon return {n,x} is sigma-clipped and zero-averaged
c - NOTE: Upon return, the number of data points is unchanged, because
c outlying items are substituted by the average of the time
c series.
c--------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333)
c
c....................................... Iterative sigma-clipping
nclip=0
1 continue
call astat(x,n,ave,sig)
si1=asig*sig
k=0
dd=-1.0d30
do 2 i=1,n
d=dabs(x(i)-ave)
if(d.lt.dd) go to 2
dd=d
k=i
2 continue
if(dd.lt.si1) go to 3
nclip=nclip+1
x(k)=ave
go to 1
3 continue
c
c....................................... Compute zero-average time series
call astat(x,n,ave,sig)
do 5 i=1,n
x(i)=x(i)-ave
5 continue
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nclip
if(ipr.eq.1) write(*,11) nclip
if(ipr.eq.2) write(*,12) nclip
10 format('>> nclip [after *PLACE UNKNOWN*] = ',i3)
11 format('>> nclip [after *datain_tfa_four*] = ',i3)
12 format('>> nclip [after *one_rec*] = ',i3)
c
c
return
end
c
c
subroutine sigma_clipping_r(ipr,asig,n,x,xf,nclip)
c------------------------------------------------------------------------
c - Performs iterative sigma-clipping of {n,x} in *respect of {xf}*
c - Level of clipping: asig*sigma
c - Upon return {n,x} is sigma-clipped
c - N is the SAME as that of the input data, because the outlying
c data points are replaced by the corresponding value of {xf}
c------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),xf(30333),y(30333)
c
c............................... Iterative sigma-clipping based on {x-xf}
nclip=0
do 4 i=1,n
y(i)=x(i)-xf(i)
4 continue
c
1 continue
call astat(y,n,ave,sig)
si1=asig*sig
k=0
dd=-1.0d30
do 2 i=1,n
d=dabs(y(i)-ave)
if(d.lt.dd) go to 2
dd=d
k=i
2 continue
if(dd.lt.si1) go to 3
nclip=nclip+1
x(k)=xf(k)
y(k)=0.0d0
go to 1
3 continue
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nclip
if(ipr.eq.1) write(*,11) nclip
if(ipr.eq.2) write(*,12) nclip
10 format(/,'>> nclip [after *PLACE UNKNOWN*] = ',i3)
11 format(/,'>> nclip [after *datain_tfa_four*] = ',i3)
12 format(/,'>> nclip [after *one_rec*] = ',i3)
c
c
return
end
c
c
subroutine ts_ar_dat(npr,n,t,x_int,xf,xfa,xw)
c---------------------------------------------------------------------
c This one prepares the data for 'ar_test_plot.sh' to test AR
c---------------------------------------------------------------------
c
c INPUT: {t} = Timebase of {x_int}
c {x_int} = Original time series [interpolated to {t}]
c {xf} = FOUR part of the FOUR+TFA fit to {x_int}
c {xfa} = Final fit for the smooth stellar variability,
c employing AR modeling at the edges and Fourier
c elsewhere (AR modeling is a part of the
c FOUR+TFA modeling)
c {xw} = Residual of the complete model:
c {xw}={x_int}-{xfa}-{tfats}
c
c OUTPUT: 'ts_ar.dat'
c
c NOTES: - If isel=1, then, depending on the decision made by
c routine "ar_vs_four", the {xfa} array may be mixed
c (i.e., FOUR on one side and AR on the other side)
c
c=====================================================================
c
c NOTES:
c
c *Because the ZERO point of the input time series {x} is fitted
c by the TFA vector, the FOUR/AR fits could be shifted by a
c constant. For this reason, we optimally shift {xf} and {xfa}
c ONLY in this routine.
c
c *X_fit=1 indicates where the Gibbs oscillations have been
c attempted to be remedied.
c
c *We assume that the moments of time are monotonic
c (i.e., t(i) > t(i-1))
c
c=====================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x_int(30333),xf(30333),xfa(30333),xw(30333)
c
tmin=t(1)
a_x=0.0d0
a_xf=0.0d0
a_xfa=0.0d0
do 1 i=1,n
a_x=a_x+x_int(i)
a_xf=a_xf+xf(i)
a_xfa=a_xfa+xfa(i)
1 continue
a_x=a_x/dfloat(n)
a_xf=a_xf/dfloat(n)-a_x
a_xfa=a_xfa/dfloat(n)-a_x
c
open(1,file='ts_ar.dat')
write(1,44)
44 format('#',/,
&'# TIME X_input X_fit X_pred X_Four',
&' X_W',/,
&'#-----------------------------------------------------',
&'------------')
n1=npr+1
n2=n-npr
t1=t(1)
do 43 i=1,n
x_fit=1.0d0
if(i.lt.n1) go to 60
if(i.gt.n2) go to 60
x_fit=xfa(i)-a_xfa
60 continue
write(1,45) t(i)-tmin,x_int(i),x_fit,xfa(i)-a_xfa,
&xf(i)-a_xf,xw(i)
45 format(6(1x,f10.6))
43 continue
close(1)
c
return
end
c
c
subroutine end_pol(npol,npr,n,t,y,xfa)
c---------------------------------------------------------------------
c Robust fit of low-order polynomials to the AR extrapolations
c of the end sections
c=====================================================================
c NOTE: For continuity/stability [avoiding edge effect in the fit],
c we fit a larger chunk than the one given by "npr". However,
c we substitute only "npr" values with the data free from
c the polynomial trend.
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),y(30333),xfa(30333)
dimension u(30333),v(30333),xpi(30333)
c
if(npol.eq.0) return
c
nn=idint(1.5d0*dfloat(npr))
c
c............................................ RIGHT end
n1=n-nn+1
n2=n
k=0
do 80 i=n1,n2
k=k+1
u(k)=t(i)
v(k)=y(i)-xfa(i)
80 continue
call robust_pol(nn,npol,u,v,xpi,sig,stdv)
k=0
n3=n-npr+1
do 81 i=n1,n2
k=k+1
if(i.lt.n3) go to 81
xfa(i)=xfa(i)+xpi(k)
81 continue
c
c............................................ LEFT end
n1=1
n2=nn
k=0
do 1 i=n1,n2
k=k+1
u(k)=t(i)
v(k)=y(i)-xfa(i)
1 continue
call robust_pol(nn,npol,u,v,xpi,sig,stdv)
k=0
n3=npr
do 2 i=n1,n2
k=k+1
if(i.gt.n3) go to 2
xfa(i)=xfa(i)+xpi(k)
2 continue
c
return
end
c
c
subroutine robust_pol(n,np,t,x,xi,sig,stdv)
c--------------------------------------------------------------
c This is a robust POLYNOMIAL fit with Cauchy weights
c==============================================================
c
c DATE: 2020.06.03
c
c NOTE: The constant in the weight factor is given by:
c (s3*RMS)**2, where RMS is the RMS of the past fit
c--------------------------------------------------------------
c
c INPUT: n = Number of data points
c np = Polynomial order
c {t} = Time values
c {x} = Time series values
c
c OUTPUT: {xi} = {best fitting polynomial}
c sig = Standard deviation of {x}-{pf}
c stdv = Non-weighted and unbiased standard deviation
c of {x}-{pf}
c
c==============================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),xw(30333),we(30333)
dimension ti(30333),xi(30333),c(995)
c
nd = np
itp = 1
isp = 1
nn = n
c
s3=1.0d0
nit=20
sig=1.0d30
sig0=1.0d0
eps=0.001d0
call astat(x,n,aver,sig0)
c
c D = SUM (x-f)^2/(cc+(x-f)^2)
c
c................................ Initialize residuals
do 2 i=1,n
xw(i)=0.0d0
2 continue
c
c................................ START ITERATION CYCLE
do 1 it=1,nit
c
if(it.le.2) go to 4
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 1
sig0=sig
4 continue
c
c................................ Update weights
cc=s3*sig0
c2=cc*cc
s=0.0d0
do 5 i=1,n
dd=xw(i)
d2=dd*dd
d2=1.0d0/(c2+d2)
s=s+d2
we(i)=d2
5 continue
s=1.0d0/s
do 6 i=1,n
we(i)=s*we(i)
6 continue
c
call polfit_w(itp,isp,n,t,x,we,nd,nn,ti,xi,c)
s=0.0d0
a=0.0d0
do 8 i=1,n
d=x(i)-xi(i)
xw(i)=d
d2=d*d
s=s+we(i)*d2
a=a+d2
8 continue
sig=dsqrt(s)
stdv=dsqrt(a/dfloat(n-np))
c
1 continue
c
return
end
c
c
subroutine robust_pol2(np,n,t,x,c,sig)
c-----------------------------------------------------------------
c This is a robust POLYNOMIAL fit with Cauchy kernel
c=================================================================
c
c INPUT: np = Polynomial order
c n = Number of data points
c {t} = Time values
c {x} = Time series values
c
c OUTPUT: {c} = Coefficients of the best-fitting polynomial
c sig = Square root of the sum of the weighted
c squared differences between {x} and the
c polynomial fit. The sum of the weights is
c equal to 1.
c
c NOTE: The constant in the weight factor is given by:
c (s3*RMS)**2, where RMS is the RMS of the past fit
c
c=================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension xw(30333),we(30333)
dimension ti(30333),xi(30333),c(995)
c
r=2.0d0/(t(n)-t(1))
tmid=(t(1)+t(n))/2.0d0
c
nd = np
itp = 1
isp = 1
nn = n
c
s3=1.0d0
nit=20
sig=1.0d30
sig0=1.0d0
eps=0.001d0
call astat(x,n,aver,sig0)
c
c D = SUM (x-f)^2/(cc+(x-f)^2)
c
c................................ Initialize residuals
do 2 i=1,n
xw(i)=0.0d0
2 continue
c
c................................ START ITERATION CYCLE
do 1 it=1,nit
c
if(it.le.2) go to 4
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 1
sig0=sig
4 continue
c
c................................ Update weights
cc=s3*sig0
c2=cc*cc
s=0.0d0
do 5 i=1,n
dd=xw(i)
d2=dd*dd
d2=1.0d0/(c2+d2)
s=s+d2
we(i)=d2
5 continue
s=1.0d0/s
do 6 i=1,n
we(i)=s*we(i)
6 continue
c
call polfit_w(itp,isp,n,t,x,we,nd,nn,ti,xi,c)
s=0.0d0
do 8 i=1,n
d=x(i)-xi(i)
xw(i)=d
s=s+we(i)*d*d
8 continue
sig=dsqrt(s)
c
1 continue
c
return
end
c
c
subroutine polfit_w(itp,isp,n,t,x,we,nd,nn,ti,xi,c)
c
c This is a polynomial fitting routine with *weights*
c
c Input time series to be fitted:
c
c (t(i),x(i)), i=1,2,...,n
c
c c(i), i=1,2,...,nd+1 regression coefficients
c
c nd = degree of the polynomial to be fitted ( must be lower than 20 )
c
c Output time series:
c (ti(i),xi(i)), i=1,2,...,nn
c
c ti(1)=t(1)
c ti(nn)=t(n)
c
c Input parameters/arrays: n,t,x,nd,nn
c ~~~~~~~~~~~~~~~~~~~~~~~~
c
c Output parameters/arrays: tt,xx,c
c ~~~~~~~~~~~~~~~~~~~~~~~~~
c
c Note: We assume that the moments of time are monotonic, i.e.,
c ~~~~~ t(i) > t(i-1)
c
c-------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
dimension t(30333),g(995,996)
dimension x(30333),c(995),td(21,30333)
dimension ti(30333),xi(30333),u(30333)
dimension we(30333)
c
if(n.gt.30333) write(*,*) 'Too many data points!'
if(n.gt.30333) return
if(nd.gt.20) write(*,*) 'Too high degree polynomial !'
if(nd.gt.20) return
c
np=nd+1
c
c Return, if the number of data points is too low
c
if(n.le.np) write(*,*) 'Too few data points !'
if(n.le.np) return
c
it = itp
is = isp
c
c it = 0 --- use array t(i) as the independent variable
c = 1 --- transform t(i) to the interval [-1,1]
c
c is = 0 --- compute equidistant interpolated time series
c = 1 --- use the same values of t(i) for interpolation
c (ti(i)=t(i), i=1,2,...,n)
c
if(is.ne.1) go to 1
if(n.ne.nn) write(*,100)
100 format('In the polinomial fit: nn must be equal to n if is=1!')
nn=n
1 continue
c
c transform t
c
tmin=t(1)
tmax=t(n)
r=2.0d0/(tmax-tmin)
tmid=(tmin+tmax)/2.0d0
if(it.gt.0) go to 20
do 23 i=1,n
u(i)=t(i)
23 continue
go to 21
20 continue
do 22 i=1,n
u(i)=r*(t(i)-tmid)
22 continue
21 continue
c
c compute the various power of {t}
c
np1=np+1
do 500 i=1,n
td(1,i)=1.0d0
do 501 j=2,np
k=j-1
td(j,i)=u(i)*td(k,i)
501 continue
500 continue
c
c************************
c compute the fit *
c************************
c
do 1001 j=1,np
c(j)=0.0d0
1001 continue
c
c**** computation of the normal matrix
c
do 11 i=1,np
do 12 j=i,np
s=0.0d0
do 13 k=1,n
s=s+we(k)*td(i,k)*td(j,k)
13 continue
g(i,j)=s
12 continue
11 continue
do 53 i=1,np
do 54 j=i,np
g(j,i)=g(i,j)
54 continue
53 continue
c
do 55 i=1,np
s=0.0d0
do 14 k=1,n
s=s+we(k)*td(i,k)*x(k)
14 continue
g(i,np1)=s
55 continue
c
c**** compute the solution
c
call gelim1(g,c,np,IFLAG)
c
c compute the fitted/interpolated time series
c
if(is.eq.1) go to 24
dt=(tmax-tmin)/dfloat(nn-1)
do 26 i=1,nn
ti(i)=tmin+dt*dfloat(i-1)
26 continue
go to 27
24 continue
do 25 i=1,n
ti(i)=t(i)
25 continue
27 continue
c
if(it.eq.0) go to 28
do 29 i=1,nn
u(i)=r*(ti(i)-tmid)
29 continue
go to 30
28 continue
do 31 i=1,nn
u(i)=ti(i)
31 continue
30 continue
c
do 520 i=1,nn
time=u(i)
s=c(nd+1)
do 521 j=1,nd
j1=nd+1-j
s=s*time+c(j1)
521 continue
xi(i)=s
520 continue
c
return
end
c
c
subroutine cut_ts(npr,nt,t,x,n,u,v)
c----------------------------------------------------------------
c This routine cuts *npr* data points on both sides of {t,xt}
c----------------------------------------------------------------
c
c INPUT: - npr :: Number of data points to be cut
c - nt :: Original number of data points
c - {t,x} :: Original time series
c
c OUTPUT: - n :: New number of data points
c - {u,v} :: Truncated time series
c================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension u(30333),v(30333)
c
n1=npr+1
n2=nt-npr
k=0
do 4 i=1,nt
if(i.lt.n1) go to 4
if(i.gt.n2) go to 4
k=k+1
u(k)=t(i)
v(k)=x(i)
4 continue
n=k
c
return
end
c
c
subroutine AR_LR(m,n,x,c)
c------------------------------------------------------------------
c This routine computes TWO-SIDE AUTOREGRESSIVE fit
c------------------------------------------------------------------
c
c DATE: - 2020.02.17
c
c INPUT: - m :: Order of the AR model
c - {x(i); i=1,2,...,n} :: Input time series
c
c OUTPUT: - {c(j); j=1,2,...,m}
c
c NOTE: - This model is two-sided, i.e.:
c x(i) = 0.5*(c1*x(i+1)+c2*x(i+2)+ ...
c c1*x(i-1)+c2*x(i-2)+ ...
c==================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333)
dimension c(995)
dimension h(995,30333)
c
c------------------------------------------------------------------
c
c............................... Set LEFT only fit functions {h}
n1=1
n2=m
do 1 k=n1,n2
do 2 i=1,m
i2=k+i
h(i,k)=x(i2)
2 continue
1 continue
c
c............................... Set RIGHT only fit functions {h}
n1=n-m+1
n2=n
do 3 k=n1,n2
do 4 i=1,m
i1=k-i
h(i,k)=x(i1)
4 continue
3 continue
c
c............................... Set TWO-SIDED fit functions {h}
n1=m+1
n2=n-m
do 5 k=n1,n2
do 6 i=1,m
i1=k-i
i2=k+i
h(i,k)=0.5d0*(x(i1)+x(i2))
6 continue
5 continue
c
c............................... Solve the LS problem
n1=1
n2=n
call solve_norm(m,n1,n2,h,x,c)
c
return
end
c
c
subroutine LR_pred(m,npr,n,t,x,xf1,xf2,xp1,xp2,
&tp1,tp2,xw,sig)
c---------------------------------------------------------------------
c This routine extends the time series {n,t,x} in BOTH directions
c---------------------------------------------------------------------
c
c INPUT: - m :: AR order
c - npr :: Number of data points to be predicted
c - n :: Number of data points of the input time series
c - {t} :: Moments of time of the input time series
c - {x} :: Values of the input time series
c
c OUTPUT:- {xf1} :: Backward prediction in [1,n-m]
c - {xf2} :: Forward prediction in [m+1,n]
c - {xp1} :: Backward prediction: npr values, BEFORE t(1)
c - {tp1} :: Time values for {xp1}
c - {xp2} :: Forward prediction: npr values, AFTER t(n)
c - {tp2} :: Time values for {xp2}
c - {xw} :: Residuals for {x(i); i=1,2,...,n}
c - sig :: Standard deviation of the FORWARD fit
c
c
c NOTES: - ORIGINAL plan:
c The extension is made ITERATIVELY within the
c framework of Fahlman & Ulrich (1982, MN, 199, 53):
c this means prediction from the AR model and LS
c computation of the predicted values by using the
c AR coefficients computed from the observed data.
c
c - ACTUAL implementation:
c A SIMPLE AR prediction. [For one-sided prediction the
c method of Fahlman & Ulrich may be inefficient]
c
c - The input data are interpolated to an equidistant grid.
c The prediction is performed on this grid. The predicted
c values are given on this EQUIDISTANT timebase.
c
c - We assume that the moments of time are monotonic, i.e.,
c t(i) > t(i-1)
c
c=====================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension xd(30333),yd(30333)
dimension xp1(30333),xp2(30333)
dimension tp1(30333),tp2(30333)
dimension xf1(30333),xf2(30333)
dimension t(30333),x(30333)
dimension u(30333),v(30333),xw(30333)
dimension c(995)
c
write(*,8)
8 format('... AR prediction')
c
c.................................... Compute average sampling rate
tmin=t(1)
tmax=t(n)
dt=(tmax-tmin)/dfloat(n-1)
c
c.................................... Interpolate data on that grid
IQ=1
N1=n
N2=n
do 5 i=1,n
u(i)=t(i)
v(i)=x(i)
xd(i)=tmin+dt*dfloat(i-1)
5 continue
CALL INTER2(IQ,u,v,N1,xd,yd,N2)
c
c.................................... Compute AR coefficients
call AR_LR(m,n,yd,c)
c
c.................................... Backward (L) prediction from the above model
ip=-1
call ar_pred(ip,m,n,yd,c,npr,xp1)
do 6 i=1,npr
tp1(i)=tmin-dt*dfloat(i)
6 continue
c
c.................................... Forward (R) prediction from the above model
ip=1
call ar_pred(ip,m,n,yd,c,npr,xp2)
do 7 i=1,npr
tp2(i)=tmax+dt*dfloat(i)
7 continue
c
c.................................... Backward (L) data prediction
n1=1
n2=n-m
do 1 i=n1,n2
s=0.0d0
do 2 j=1,m
s=s+c(j)*yd(j+i)
2 continue
xf1(i)=s
xw(i)=s-yd(i)
1 continue
c
c.................................... Forward (R) data prediction
n1=m+1
n2=n
sum=0.0d0
do 3 i=n1,n2
s=0.0d0
do 4 j=1,m
s=s+c(j)*yd(i-j)
4 continue
xf2(i)=s
d=s-yd(i)
xw(i)=d
sum=sum+d*d
3 continue
sig=dsqrt(sum/dfloat(n2-n1+1-m))
c
return
end
c
c
subroutine ar_pred(ip,m,n,x,c,npr,xp)
c------------------------------------------------------------------
c This is an AR prediction of "npr" values preceeding (ip=-1)
c or following (ip=+1) the dataset in [1,n]
c------------------------------------------------------------------
c
c NOTE: - The index of {xp} grows as we depart from the
c 1st/last data point of the original time series
c
c==================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),xp(30333)
dimension c(995),y(995),z(995)
c
c............................... Set prediction basis {y}
if(ip.eq.-1) go to 6
do 3 j=1,m
k=n-j+1
y(j)=x(k)
3 continue
go to 7
6 continue
do 8 j=1,m
y(j)=x(j)
8 continue
7 continue
c
c............................... Start predicting "npr" values
do 1 i=1,npr
c
s=0.0d0
do 2 j=1,m
s=s+c(j)*y(j)
2 continue
xp(i)=s
c
c............................... Copy shifted {y} to {z}
do 5 j=2,m
z(j)=y(j-1)
5 continue
z(1)=s
c............................... Update {y}
do 4 j=1,m
y(j)=z(j)
4 continue
c
1 continue
c
return
end
c
c
subroutine solve_norm(m,n1,n2,h,x,c)
c------------------------------------------------------------------
c This routine solves the Least Squares problem derived from
c the design matrix {h(i,j)} and data array {x(k)}:
c
c x_est (k) = SUM_j c(j)*h(j,k)
c------------------------------------------------------------------
c
c DATE: - 2020.01.31
c
c INPUT: - m :: Number of parameters to be fitted
c - n1 :: Index of the first array element
c to be fitted
c - n2 :: Index of the last array element
c to be fitted
c - {h(i,j)} :: m x (n2-n1+1) matrix
c - {x} :: Data to be fitted
c
c OUTPUT: - {c} :: Regression coefficients
c
c NOTE: - The design matrix {h} is defined in [n1,n2], where
c n2 > n1 and 1 < n1,n2 < n
c
c==================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),c(995)
dimension h(995,30333)
dimension g(995,996)
c
c------------------------------------------------------------------
c
m1=m+1
c
c............................... Compute normal matrix
do 3 i=1,m
do 4 j=i,m
s=0.0d0
do 5 k=n1,n2
s=s+h(i,k)*h(j,k)
5 continue
g(i,j)=s
g(j,i)=s
4 continue
3 continue
c
c............................... Compute RHS
do 7 i=1,m
s=0.0d0
do 15 k=n1,n2
s=s+h(i,k)*x(k)
15 continue
g(i,m1)=s
7 continue
c
c............................... Compute solution
call gelim1(g,c,m,IFLAG)
c
c
return
end
c
c
subroutine ar_four_ar(npr,nt,t,xf,tp1,xp1,tp2,xp2,xfa)
c----------------------------------------------------------------------
c This routine combines the FOURIER and the AR time series
c----------------------------------------------------------------------
c
c INPUT: - npr :: Number of AR-predicted data points
c - nt :: Original number of data points
c - {t} :: Moments of time of the original time series
c - {xf} :: Values of the FOURIER fit on {t}
c - {tp1} :: Time values for {xp1}
c - {xp1} :: Backward prediction: BEFORE t(npr+1)
c - {tp2} :: Time values for {xp2}
c - {xp2} :: Forward prediction: AFTER t(nt-npr)
c
c OUTPUT: - {xfa} :: {xp1}+{xf}+{xp2} combined time series on {t}
c
c NOTES: - {xf} is kept ONLY in the "inner part" of the time
c series (i.e., between i=npr+1 and nt-npr
c - The predictions are given on EQUIDISTANT timebase.
c At the end, these values are INTERPOLATED back to the
c INPUT timebase given by {t}.
c
c======================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),xf(30333),xfa(30333)
dimension xp1(30333),xp2(30333)
dimension tp1(30333),tp2(30333)
dimension u(30333),v(30333)
dimension xd(30333),yd(30333)
c
c....................................... INTERPOLATE prediction
c back to timebase {t}
IQ=1
N1=npr
N2=npr
do 1 i=1,npr
j=npr-i+1
u(i)=tp1(j)
v(i)=xp1(j)
xd(i)=t(i)
1 continue
CALL INTER2(IQ,u,v,N1,xd,yd,N2)
do 2 i=1,npr
xfa(i)=yd(i)
2 continue
do 3 i=1,npr
u(i)=tp2(i)
v(i)=xp2(i)
xd(i)=t(i+nt-npr)
3 continue
CALL INTER2(IQ,u,v,N1,xd,yd,N2)
do 4 i=1,npr
xfa(i+nt-npr)=yd(i)
4 continue
c
c............................... Fill in the non-extrapolated part
c by the FOURIER fit
n1=npr+1
n2=nt-npr
do 5 i=n1,n2
xfa(i)=xf(i)
5 continue
c
return
end
c
c
subroutine save_rec(n,t,x_rec,xft,xsyn,trsyn,
&tfa_r,four_r,qtran)
c----------------------------------------------------------
c Save reconstructed/filtered time series {x_rec} with
c the best trapezoidal fit {xft} and the synthetic
c signal {xsyn}, if test data are used [if not, then
c this array is set equal to zero]
c ... Also saved are: {tfa_r} and {four_r}, TFA and
c FOUR vectors from the full/reconstructed model.
c----------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x_rec(30333)
dimension xft(30333),xsyn(30333),trsyn(30333)
dimension tfa_r(30333),four_r(30333)
c
open(18,file='rec_tfa.dat')
write(18,56) qtran
56 format('#',/,'# QTRAN = ',f7.5,/,
&'#',91('='),/,
&'# TIME',11x,'X_REC',9x,'XFT',10x,'SYNT',7x,'SYN_TRAN',
&' FOUR_REC TFA_REC',/,
&'#',91('-'))
do 1 i=1,n
write(18,2) t(i),x_rec(i),xft(i),xsyn(i),trsyn(i),
&four_r(i),tfa_r(i)
1 continue
2 format(f14.7,6(2x,f11.7))
close(18)
c
return
end
c
c
subroutine save_sp(nb_sp,nf1,fmin1,df1,p)
c----------------------------------------------------------
c - This routine saves the "nb_sp-binned-maximum-value"
c spectrum as "sp.dat"
c----------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension p(155555)
c
nst = nf1/nb_sp
nf2 = nf1/nst
df2 = nst*df1
c
open(25,file='sp.dat')
do 1 i=1,nf1,nst
i1=i
f1=fmin1+df1*dfloat(i1-1)
i2=i+nst
f2=fmin1+df1*dfloat(i2-1)
freq=(f1+f2)*0.5d0
s=-1.0d30
i0=i1-1
do 2 j=1,nst
k=i0+j
s=dmax1(s,p(k))
2 continue
write(25,3) freq,s
3 format(f10.6,1x,f12.8)
1 continue
close(25)
c
return
end
c
c
subroutine save_spec(sname,nb_sp,nf1,fmin1,df1,p)
c----------------------------------------------------------
c - This routine saves the "nb_sp-binned-maximum-value"
c spectrum as "spec/sname.sp"
c----------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 sname
character*17 spname
c
dimension p(155555)
c
spname= trim('spec/'//sname//'.sp')
c
nst = nf1/nb_sp
nf2 = nf1/nst
df2 = nst*df1
c
open(19,file=spname)
do 1 i=1,nf1,nst
i1=i
f1=fmin1+df1*dfloat(i1-1)
i2=i+nst
f2=fmin1+df1*dfloat(i2-1)
freq=(f1+f2)*0.5d0
s=-1.0d30
i0=i1-1
do 2 j=1,nst
k=i0+j
s=dmax1(s,p(k))
2 continue
write(19,3) freq,s
3 format(f10.6,1x,f12.8)
1 continue
close(19)
c
return
end
c
c
subroutine save_sp_iw(nb_sp,nf1,fmin1,df1,p,spn)
c
c----------------------------------------------------------
c This routine saves the "nbsp-binned-maximum-value"
c frequency spectrum in the working directory under
c the name "spn"
c---------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*8 spn
c
dimension p(155555)
c
nst = nf1/nb_sp
nf2 = nf1/nst
df2 = nst*df1
c
open(19,file=spn)
do 1 i=1,nf1,nst
i1=i
f1=fmin1+df1*dfloat(i1-1)
i2=i+nst
f2=fmin1+df1*dfloat(i2-1)
freq=(f1+f2)*0.5d0
s=-1.0d30
i0=i1-1
do 2 j=1,nst
k=i0+j
s=dmax1(s,p(k))
2 continue
write(19,3) freq,s
3 format(f10.6,1x,f12.8)
1 continue
close(19)
c
return
end
c
c
subroutine save_lc_iw(n,t,y,xbt,lcn)
c
c----------------------------------------------------------
c This routine saves the time series {n,t,y} in the
c working directory under the name "lcn"
c---------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*8 lcn
c
dimension t(30333),y(30333),xbt(30333)
c
open(19,file=lcn)
do 1 i=1,n
write(19,3) t(i),y(i),xbt(i)
3 format(f13.7,2(1x,f10.7))
1 continue
close(19)
c
return
end
c
c
subroutine save_pa_iw(sname,p0,epc,pan)
c-------------------------------------------------------------------
c This routine saves folding parameters
c-------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*8 pan
character*9 sname
c
open(24,file=pan)
nper=1
f0=1.0d0/p0
write(24,40) sname,nper,f0,epc
40 format(a9,/,i2,/,f11.7,/,f13.6)
close(24)
c
return
end
c
c
subroutine parin(isurvey,jsm,itemp,nstar,fmin1,fmax1,nf1,
&iwhite,nb,qmi,qma,itest,tfreq0,tfreq1,dep0,dep1,t14_0,
&t14_1,pht0,pht1,frin0,frin1,amootv,phootv,frootv,iflare,
&igen,std,fdep,iseed,issp,isflc,ispl,iplot,asig,ic_flare,
&ibg,ixy,m_four,m_opt)
c-------------------------------------------------------------------
c This routine reads in parameters more frequently changed while
c running this code
c===================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 name
character*1 eq
c
call SYSTEM('grep -v "#" tran_k2_v0.par > temp.dat')
open(1,file='temp.dat')
c
read(1,*) name,eq,isurvey
read(1,*) name,eq,jsm
read(1,*) name,eq,itemp
read(1,*) name,eq,nstar
read(1,*) name,eq,fmin1
read(1,*) name,eq,fmax1
read(1,*) name,eq,nf1
read(1,*) name,eq,iwhite
read(1,*) name,eq,nb
read(1,*) name,eq,qmi
read(1,*) name,eq,qma
c=================================
read(1,*) name,eq,itest
read(1,*) name,eq,tfreq0
read(1,*) name,eq,dep0
read(1,*) name,eq,t14_0
read(1,*) name,eq,pht0
read(1,*) name,eq,frin0
c=================================
read(1,*) name,eq,tfreq1
read(1,*) name,eq,dep1
read(1,*) name,eq,t14_1
read(1,*) name,eq,pht1
read(1,*) name,eq,frin1
c=================================
read(1,*) name,eq,frootv
read(1,*) name,eq,amootv
read(1,*) name,eq,phootv
read(1,*) name,eq,iflare
read(1,*) name,eq,igen
read(1,*) name,eq,std
read(1,*) name,eq,fdep
read(1,*) name,eq,iseed
c=================================
read(1,*) name,eq,issp
read(1,*) name,eq,isflc
read(1,*) name,eq,ispl
read(1,*) name,eq,iplot
read(1,*) name,eq,asig
read(1,*) name,eq,ic_flare
read(1,*) name,eq,ibg
read(1,*) name,eq,ixy
read(1,*) name,eq,m_four
read(1,*) name,eq,m_opt
c=================================
close(1)
c
return
end
c
c
subroutine minmax(n,x,xmin,xmax)
c
c This routine computes (min,max) of {x(i); i=1,...,n}
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333)
c
xmin = 1.0d30
xmax =-1.0d30
do 1 i=1,n
xx=x(i)
xmin = dmin1(xmin,xx)
xmax = dmax1(xmax,xx)
1 continue
c
return
end
c
c
FUNCTION ARAND(IX)
implicit real*8 (a-h,o-z)
INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K
DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/
XHI = IX/B16
XALO = (IX - XHI * B16) * A
LEFTLO = XALO / B16
FHI = XHI * A + LEFTLO
K = FHI / B15
IX = ((( XALO - LEFTLO * B16 ) - P ) + (FHI - K * B15 ) * B16) + K
IF (IX.LT.0) IX = IX + P
ARAND = DFLOAT(IX) * 4.65661287D-10
RETURN
END
C
C
FUNCTION GAUSS(ISEED)
C
C this routine generates Gaussian random variables of unit variance
C
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
C
C******** TRANSFORMING TO GAUSSIAN VARIABLE *******
C******** *******
C SOURCE: NUMERICAL RECIPES, P. 203
C
1 CONTINUE
P1=2.0D0*ARAND(ISEED)-1.0D0
P2=2.0D0*ARAND(ISEED)-1.0D0
RR=P1*P1+P2*P2
IF(RR.GE.1.0D0) GO TO 1
FAC=DSQRT(-2.0D0*DLOG(RR)/RR)
GAUSS=P1*FAC
C
RETURN
END
c
c
subroutine astat(y,nn,aver,disper)
c
c This routine computes the average and unbiased standard deviation of {y}
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
dimension y(30333)
c
s=0.0d0
do 1 i=1,nn
s=s+y(i)
1 continue
aver=s/dfloat(nn)
di=0.0d0
do 2 i=1,nn
d=y(i)-aver
di=di+d*d
2 continue
disper=0.0d0
if(nn.gt.1) disper=dsqrt(di/dfloat(nn-1))
c
return
end
C
C
subroutine mfourdec(epoch,n,t,x,m,w,xw,sig,a00,aa,pp)
C
C----------------------------------------------------------------
C THIS CODE FITS A FOURIER SUM TO A TIME-SERIES
C----------------------------------------------------------------
C
C EPOCH = THE FIT IS COMPUTED ON A TIMEBASE OF {T(I)-EPOCH}
C N = NUMBER OF DATA POINTS
C T(I) = TIME
C X(I) = TIME-SERIES
C M = ORDER OF THE FOURIER FIT
c W(J) = J-TH FREQUENCY (J=1,2,...,M)
C XW(I) = PREWHITENED TIME-SERIES
C SIG = STANDARD DEVIATION OF THE FIT [UNBAISED ESTIMATE]
C A00 = ZERO FREQUENCY COMPONENT
C AA(J) = AMPLITUDE OF THE J-TH COMPONENT
C PP(J) = PHASE OF THE J-TH COMPONENT (IN [RAD])
C
C FOURIER DECOMPOSITION HAS THE FOLLOWING FORM:
C
C A0 + A1*SIN(2*PI*F1*(T-EPOCH)+FI1) + ...
C
C================================================================
C
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER*4 (I-N)
C
DIMENSION X(30333),T(30333),G(995,996),W(995)
DIMENSION U(30333),XW(30333)
DIMENSION A(995),AA(995),PP(995)
C
P7 = 8.0D0*DATAN(1.0D0)
c
DO 30 I=1,N
U(I)=T(I)-EPOCH
30 CONTINUE
C
M2=2*M
M3=M2+1
M4=M3+1
CALL NORM(N,M,U,X,W,G)
CALL gelim1(G,A,M3,IFLAG)
A0=A(1)
C
C Compute amplitudes and phases
C
K=0
DO 37 J=2,M3,2
K=K+1
J1=J+1
AMP=DSQRT(A(J)*A(J)+A(J1)*A(J1))
PH=A(J1)/AMP
PH=DACOS(PH)
IF(A(J).LT.0.0D+00) PH=P7-PH
A(K)=AMP
AA(K)=AMP
PP(K)=PH
37 CONTINUE
A00=A0
C
C Compute residuals and standard deviation
C
DISP=0.0d0
DO 2000 I=1,N
S=A00
TIME=U(I)
DO 2001 K=1,M
PHASE=W(K)*TIME
S=S+AA(K)*DSIN(P7*PHASE+PP(K))
2001 CONTINUE
DEVI=X(I)-S
XW(I)=DEVI
DISP=DISP+DEVI*DEVI
2000 CONTINUE
SIG=DSQRT(DISP/DFLOAT(N-2*M-1))
C
RETURN
END
C
C
SUBROUTINE NORM(N,M,T,X,F,G)
C CALCULATION OF THE ELEMENTS OF THE NORMAL MATRIX AND R.H.S.
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
DIMENSION G(995,996),T(30333),X(30333),F(995)
P7=8.0D0*DATAN(1.0D0)
M3=2*M+1
M4=M3+1
L=0
G(1,1)=DFLOAT(N)
DO 1 J=2,M3,2
S=0.0D+00
C=0.0D+00
L=L+1
J1=J+1
FL=F(L)
DO 2 I=1,N
PH=T(I)*FL
PH=P7*PH
S=S+DSIN(PH)
2 C=C+DCOS(PH)
G(1,J)=C
1 G(1,J1)=S
L1=0
DO 3 J=2,M3,2
L1=L1+1
J1=J+1
L2=L1-1
FL1=F(L1)
DO 4 K=J,M3,2
K1=K+1
L2=L2+1
FL2=F(L2)
S1=0.0D+00
S2=0.0D+00
C1=0.0D+00
C2=0.0D+00
DO 5 I=1,N
PH1=T(I)*FL1
PH1=P7*PH1
PH2=T(I)*FL2
PH2=P7*PH2
SPH1=DSIN(PH1)
CPH1=DCOS(PH1)
SPH2=DSIN(PH2)
CPH2=DCOS(PH2)
C1=C1+CPH1*CPH2
S1=S1+CPH2*SPH1
C2=C2+SPH2*CPH1
5 S2=S2+SPH1*SPH2
G(J,K)=C1
G(J,K1)=C2
G(J1,K)=S1
4 G(J1,K1)=S2
3 CONTINUE
S=0.0D+00
DO 6 I=1,N
6 S=S+X(I)
G(1,M4)=S
L=0
DO 7 J=2,M3,2
S=0.0D+00
C=0.0D+00
L=L+1
J1=J+1
FL=F(L)
DO 8 I=1,N
PH=T(I)*FL
PH=P7*PH
S=S+X(I)*DSIN(PH)
8 C=C+X(I)*DCOS(PH)
G(J,M4)=C
7 G(J1,M4)=S
DO 9 J=1,M3
DO 10 K=J,M3
10 G(K,J)=G(J,K)
9 CONTINUE
RETURN
END
C
c
subroutine fdft(n,t,y,nf,fmin,df,p,pfr,pam)
c
c-------------------------------------------------------------
c This is a FAST single frequency DFT routine based on
c the idea of Kurtz (1985, MN, 213, 773)
c [and mine ... a bit earlier ... unpublished ...]
c-------------------------------------------------------------
c
c {t(i),y(i)}, i=1,2, ... , n --- time-series
c {p(j)}, j=1,2, ... , nf --- ampl. spectr.
c fmin --- min. test frequency
c df --- frequency step
c pfr --- peak frequency
c pam --- peak amplitude
c
c WARNING: - Maximum number of data points MUST be
c ~~~~~~~~ lower than *30333*
c
c - Maximum number of spectrum points MUST be
c lower than *155555*
c.............................................................
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
real*4 u(30333),v(30333),pi2,fp,sum
real*4 ui,vi,a,da,omin,dom,s1i,c1i,s0i,c0i,s2i,sj,cj
real*4 s(155555),c(155555)
real*4 s1(30333),c1(30333),s0(30333),c0(30333)
c
dimension t(*),y(*),p(*)
c
fp = 2.0/float(n)
pi2 = 8.0*atan(1.0)
pam = -1.0
c
if(n.gt.30333) write(*,*) ' n MUST be lower that 30333!'
if(n.gt.30333) return
if(nf.gt.155555) write(*,*) ' nf MUST be lower that 155555!'
if(nf.gt.155555) return
c
c=================================
c Set temporal time series
c=================================
c
sum=0.0
t1=t(1)
do 103 i=1,n
u(i)=sngl(t(i)-t1)
sum=sngl(sum+y(i))
103 continue
sum=sngl(sum/dfloat(n))
do 109 i=1,n
v(i)=sngl(y(i)-sum)
109 continue
c
omin = sngl(pi2*fmin)
dom = sngl(pi2*df)
do 3 i=1,nf
s(i) = 0.0
c(i) = 0.0
3 continue
c
do 5 i=1,n
ui = u(i)
vi = v(i)
a = omin*ui
da = dom*ui
s1(i) = vi*sin(a)
c1(i) = vi*cos(a)
s0(i) = sin(da)
c0(i) = cos(da)
5 continue
c
c... Compute amplitude spectrum
c
do 1 i=1,n
s1i = s1(i)
c1i = c1(i)
s0i = s0(i)
c0i = c0(i)
do 2 j=1,nf
s(j) = s(j) + s1i
c(j) = c(j) + c1i
s2i = s1i
s1i = s1i*c0i+c1i*s0i
c1i = c1i*c0i-s2i*s0i
2 continue
1 continue
c
do 4 j=1,nf
sj = s(j)
cj = c(j)
a = fp*sqrt(sj*sj+cj*cj)
p(j) = a
if(a.lt.pam) go to 4
pam=a
pfr=fmin+df*float(j-1)
4 continue
c
return
end
c
c
subroutine eebls(n,t,x,u,v,nf,fmin,df,nb,qmi,qma,
&p,bper,bpow,depth,qtran,in1,in2)
c
c------------------------------------------------------------------------
c >>>>>>>>>>>> This routine computes BLS spectrum <<<<<<<<<<<<<<
c
c [ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ]
c
c This is the slightly modified version of the original BLS routine
c by considering Edge Effect (EE) as suggested by
c Peter R. McCullough [ pmcc@stsci.edu ].
c
c This modification was motivated by considering the cases when
c the low state (the transit event) happened to be devided between
c the first and last bins. In these rare cases the original BLS
c yields lower detection efficiency because of the lower number of
c data points in the bin(s) covering the low state (the transit event).
c
c For further comments/tests see www.konkoly.hu/staff/kovacs.html
c
c The following additional modifications were made on June 7, 2004:
c
c - Constraint concerning the minimum number of bins in the transit
c has been abandoned.
c - Total time span of the time series is computed without the
c assumption of sorted time values.
c - Relative transit length is computed by using the bin indices of
c the ingress and egress and not as the ratio of the number of
c data points in the two states.
c------------------------------------------------------------------------
c
c Input parameters:
c ~~~~~~~~~~~~~~~~~
c
c n = number of data points
c t = array {t(i)}, containing the time values of the time series
c x = array {x(i)}, containing the data values of the time series
c u = temporal/work/dummy array, must be dimensioned in the
c calling program in the same way as {t(i)}
c v = the same as {u(i)}
c nf = number of frequency points in which the spectrum is computed
c fmin = minimum frequency (MUST be > 0)
c df = frequency step
c nb = number of bins in the folded time series at any test period
c qmi = minimum fractional transit length to be tested
c qma = maximum fractional transit length to be tested
c
c Output parameters:
c ~~~~~~~~~~~~~~~~~~
c
c p = array {p(i)}, containing the values of the BLS spectrum
c at the i-th frequency value -- the frequency values are
c computed as f = fmin + (i-1)*df
c bper = period at the highest peak in the frequency spectrum
c bpow = value of {p(i)} at the highest peak
c depth= depth of the transit at *bper*
c qtran= fractional transit length [ T_transit/bper ]
c in1 = bin index at the start of the transit [ 0 < in1 < nb+1 ]
c in2 = bin index at the end of the transit [ 0 < in2 < nb+1 ]
c
c Remarks:
c ~~~~~~~~
c
c -- *fmin* MUST be greater than *1/total time span*
c -- *nb* MUST be lower than *nbmax*.
c -- Dimensions of arrays {y(i)} and {ibi(i)} MUST be greater than
c or equal to *nbmax*.
c -- The lowest number of points allowed in the transit is equal
c to MAX(minbin,qmi*N), where *qmi* is the minimum transit
c length/trial period, *N* is the total number of data points,
c *minbin* is the preset minimum number of the data points in
c the transit.
c
c========================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(*),x(*),u(*),v(*),p(*)
dimension y(30333),ibi(30333)
c
minbin = 10
nbmax = 30333
jn1 = 1
jn2 = 1
rn3 = 0.0
s3 = 0.0
nb0 = nb
c
if(nb.gt.nbmax) write(*,*) ' NB > NBMAX !!'
if(nb.gt.nbmax) nb=nbmax
c
c... Compute total time span
c
tmin = 1.0d30
tmax = -1.0d30
do 10 i=1,n
tmin=dmin1(t(i),tmin)
tmax=dmax1(t(i),tmax)
10 continue
tot=tmax-tmin
c
fminold=fmin
if(fmin.lt.1.0d0/tot) write(*,*) ' fmin < 1/T !!'
if(fmin.lt.1.0d0/tot) fmin=1.1/tot
c------------------------------------------------------------------------
c
rn=dfloat(n)
rnb=dfloat(nb)
kma=idint(rnb*qma)+1
kkmi=idint(rn*qmi)
if(kkmi.lt.minbin) kkmi=minbin
bpow=0.0d0
c
c The following variables are defined for the extension
c of arrays ibi() and y() [ see below ]
c
nb1 = nb+1
nbkma = nb+kma
c
c=================================
c Set temporal time series
c=================================
c
s=0.0d0
do 103 i=1,n
u(i)=t(i)-tmin
s=s+x(i)
103 continue
s=s/rn
do 109 i=1,n
v(i)=x(i)-s
109 continue
c
c******************************
c Start period search *
c******************************
c
do 100 jf=1,nf
f0=fmin+df*dfloat(jf-1)
p0=1.0d0/f0
c
c======================================================
c Compute folded time series with *p0* period
c======================================================
c
do 101 j=1,nb
y(j) = 0.0d0
ibi(j) = 0
101 continue
c
do 102 i=1,n
ph = u(i)*f0
ph = ph-idint(ph)
j = 1 + idint(nb*ph)
ibi(j) = ibi(j) + 1
y(j) = y(j) + v(i)
102 continue
c
c-----------------------------------------------
c Extend the arrays ibi() and y() beyond
c nb by wrapping
c
do 104 j=nb1,nbkma
jnb = j-nb
ibi(j) = ibi(jnb)
y(j) = y(jnb)
104 continue
c-----------------------------------------------
c
c===============================================
c Compute BLS statistics for this period
c===============================================
c
power=0.0d0
c
do 1 i=1,nb
s = 0.0d0
kk = 0
nb2 = i+kma
do 2 j=i,nb2
kk = kk+ibi(j)
s = s+y(j)
if(kk.lt.kkmi) go to 2
rn1 = dfloat(kk)
pow = s*s/(rn1*(rn-rn1))
if(pow.lt.power) go to 2
power = pow
jn1 = i
jn2 = j
rn3 = rn1
s3 = s
2 continue
1 continue
c
power = dsqrt(power)
p(jf) = power
c
if(power.lt.bpow) go to 100
bpow = power
in1 = jn1
in2 = jn2
qtran = dfloat(in2-in1+1)/rnb
depth = -s3*rn/(rn3*(rn-rn3))
bper = p0
c
100 continue
c
c Edge correction of transit end index
c
if(in2.gt.nb) in2 = in2-nb
c
c Restore original value of 'fmin' and 'nb'
c
fmin=fminold
nb=nb0
c
return
end
c
c
subroutine weebls_new(n,t,x,w,u,v,nf,fmin,df,nb,qmi,qma,
&p,bper,bpow,depth,qtran,in1,in2)
c
c------------------------------------------------------------------------
c This routine computes BLS spectrum [ see A&A, Vol. 391, 369 ]
c
c Version: November 06, 2018
c ~~~~~~~~
c
c Description, modifications:
c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
c
c - Fits boxes with the aid of Least Squares principle to the folded
c time series for *nf* trial periods.
c
c - Considers Edge Effect:
c
c This was suggested by Peter R. McCullough [ pmcc@stsci.edu ].
c The modification was motivated by considering the case when
c the low state (the transit event) is accidently devided between
c the first and last bins. In these rare cases the original BLS
c yields lower detection efficiency because of the lower number
c of data points in the bin(s) covering the low state (the transit
c event). For further comments/tests see
c www.konkoly.hu/staff/kovacs.html
c
c - Uses weights in order to consider observational errors:
c
c Implementation was suggested by Kailash Sahu [ ksahu@stsci.edu ].
c Weighting enters in the general formulation of BLS but it was
c considered originally as time consuming and usually not necessary.
c
c - For more technical details and modifications relative to the
c very first version of 2002, see *Remarks* below.
c------------------------------------------------------------------------
c
c Input parameters:
c ~~~~~~~~~~~~~~~~~
c
c n = number of data points
c t = array {t(i)}, contains the time values of the time series
c x = array {x(i)}, contains the data values of the time series
c w = array {w(i)}, contains the weights for each item of data
c u = temporal/work/dummy array, must be dimensioned in the
c calling program in the same way as {t(i)}
c v = the same as {u(i)}
c nf = number of frequency points in which the spectrum is computed
c fmin = minimum frequency (MUST be > 0)
c df = frequency step
c nb = number of bins in the folded time series at any test period
c qmi = minimum fractional transit length to be tested
c qma = maximum fractional transit length to be tested
c
c Output parameters:
c ~~~~~~~~~~~~~~~~~~
c
c p = array {p(i)}, contains the values of the BLS spectrum.
c The frequency corresponding to p(i) is equal to
c f = fmin + (i-1)*df
c bper = period at the highest peak in the frequency spectrum
c bpow = value of {p(i)} at the highest peak
c depth= depth of the transit at *bper*
c qtran= fractional transit length [ T_transit/bper ]
c in1 = bin index at the start of the transit [ 0 < in1 < nb+1 ]
c in2 = bin index at the end of the transit [ 0 < in2 < nb+1 ]
c
c
c Remarks:
c ~~~~~~~~
c
c -- *fmin* MUST be greater than *1/total time span*
c -- *nb* MUST be lower than *nbmax*
c -- Dimensions of arrays {y(i)} and {wib(i)} MUST be set equal
c to 2*nbmax for complete security against array index overflow
c -- The sum of the weights in the transit cannot be smaller
c than *qmi*. This helps avoiding statistically unstable
c situations.
c -- In the original BLS code a check was made for the number of
c bins covering the transit. This is considered to be unimportant,
c and therefore omitted, because the function of this check is
c already fulfilled by a similar check of the sum of the weights
c in the transit.
c -- The fractional transit length *qtran* is computed from the
c bin indices of the ingress and egress and not from the number
c of data points in the transit, as in the earlier versions of
c BLS.
c -- The time values {t(i)} need not to be sorted.
c -- A check is made on the weights {w(i)} in order to verify if
c their sum is within 1+/-eps.
c -- The average of {v(i)=w(i)x(i)} is subtracted from {v(i)}, i.e.,
c the input time series {x(i)} need not to be zero-averaged.
c -- The sign definition of the output parameter *depth* follows
c astronomical magnitude definition (i.e., if depth<0, then
c it means dimming, lower brightness).
c -- Condition "if(j.eq.i) go to 2" is added to the inner loop
c on November, 6, 2018.
c========================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(*),x(*),w(*),u(*),v(*),p(*)
dimension y(30333),wib(30333)
dimension ibi(30333)
c
c------------------------------------------------------------------------
c... Parameter and {w(i)} normalization checks
c
eps = 1.0d-6
nbmax = 10000
nb0 = nb
c
c................... These are just dummy data initializations to stop
c the compiler whining on "uninitialized" variables
jn1 = 1
jn2 = 1
wkb = 0.0d0
sb = 0.0d0
c...................
c
if(nb.gt.nbmax) write(*,*) ' NB > NBMAX !!'
if(nb.gt.nbmax) nb=nbmax
s=0.0d0
do 200 i=1,n
s=s+w(i)
200 continue
d=dabs(s-1.0d0)
if(d.lt.eps) go to 10
write(*,*) ' |SUM w(i)-1.0| > eps !!'
s=1.0d0/s
do 11 i=1,n
w(i)=s*w(i)
11 continue
10 continue
c
c... Compute total time span, check *fmin*
c
tmin = 1.0d30
tmax = -1.0d30
do 201 i=1,n
tmin=dmin1(tmin,t(i))
tmax=dmax1(tmax,t(i))
201 continue
ftot=1.0d0/(tmax-tmin)
fmin5=fmin
if(fmin5.gt.ftot) go to 12
write(*,*) ' fmin < 1/T !!'
fmin5=ftot
12 continue
c------------------------------------------------------------------------
c
c... Set additional parameters
c
minbin = 5
rn = dfloat(n)
rnb = dfloat(nb)
kma = idint(qma*rnb)+1
kkmi= idint(rn*qmi)
if(kkmi.lt.minbin) kkmi=minbin
bpow= 0.0d0
c
c The following variables are defined for the extension
c of arrays {wib(i)} and {y(i)} [ see below ]
c
nb1 = nb+1
nbkma = nb+kma
c
c------------------------------------------------------------------------
c... Set temporal time series {u(i),v(i)}
c
c - t(i) --> u(i) > 0;
c - x(i), w(i) --> v(i)
c
swx=0.0d0
do 202 i=1,n
u(i)=t(i)-tmin
v(i)=w(i)*x(i)
swx=swx+v(i)
202 continue
swx=swx/rn
do 203 i=1,n
v(i)=v(i)-swx
203 continue
c
c******************************
c Start period search *
c******************************
c
do 100 jf=1,nf
f0=fmin5+df*dfloat(jf-1)
p0=1.0d0/f0
c
c======================================================
c Compute folded time series with period *p0*
c======================================================
c
do 101 j=1,nb
y(j) = 0.0d0
wib(j) = 0.0d0
ibi(j) = 0
101 continue
c
do 102 i=1,n
ph = u(i)*f0
ph = ph-idint(ph)
j = 1 + idint(nb*ph)
wib(j) = wib(j) + w(i)
ibi(j) = ibi(j) + 1
y(j) = y(j) + v(i)
102 continue
c
c-----------------------------------------------
c... Extend the arrays {wib(i)} and {y(i)}
c beyond *nb* by wrapping
c
do 104 j=nb1,nbkma
jnb = j-nb
wib(j) = wib(jnb)
ibi(j) = ibi(jnb)
y(j) = y(jnb)
104 continue
c-----------------------------------------------
c
c===============================================
c Compute BLS statistics for this period
c===============================================
c
power=0.0d0
c
do 1 i=1,nb
s = 0.0d0
wk = 0.0d0
kk = 0
nb2 = i+kma
do 2 j= i,nb2
wk = wk+wib(j)
kk = kk+ibi(j)
s = s+y(j)
c
if(kk.lt.kkmi) go to 2
if(j.eq.i) go to 2
c
pow = s*s/((1.0d0-wk)*wk)
if(pow.lt.power) go to 2
c
power = pow
jn1 = i
jn2 = j
wkb = wk
sb = s
c
2 continue
1 continue
c
power = dsqrt(power)
p(jf) = power
c
if(power.lt.bpow) go to 100
bpow = power
in1 = jn1
in2 = jn2
qtran = dfloat(in2-in1+1)/rnb
depth = -sb/((1.0d0-wkb)*wkb)
bper = p0
c
100 continue
c
c... Edge correction of transit end index
c
if(in2.gt.nb) in2 = in2-nb
c
nb=nb0
c
return
end
c
c
subroutine polfit(itp,isp,n,t,x,nd,nn,ti,xi,c)
c
c This is a polynomial fitting routine
c
c Input time series to be fitted:
c
c (t(i),x(i)), i=1,2,...,n
c
c c(i), i=1,2,...,nd+1 regression coefficients
c
c nd = degree of the polynomial to be fitted
c
c Output time series:
c (ti(i),xi(i)), i=1,2,...,nn
c
c ti(1)=t(1)
c ti(nn)=t(n)
c
c Input parameters/arrays: itp,isp,n,t,x,nd,nn
c ~~~~~~~~~~~~~~~~~~~~~~~~
c
c Output parameters/arrays: ti,xi,c
c ~~~~~~~~~~~~~~~~~~~~~~~~~
c
c Note: We assume that the moments of time are monotonic, i.e.,
c ~~~~~ t(i) > t(i-1)
c
c-------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
dimension t(30333),g(995,996)
dimension x(30333),c(995),td(21,30333)
dimension ti(30333),xi(30333),u(30333)
c
if(n.gt.30333) write(*,*) 'Too many data points!'
if(n.gt.30333) return
if(nd.gt.20) write(*,*) 'Too high degree polynomial !'
if(nd.gt.20) nd=20
c
np=nd+1
c
it = itp
is = isp
c
c it = 0 --- use array t(i) as the independent variable
c = 1 --- transform t(i) to the interval [-1,1]
c
c is = 0 --- compute equidistant interpolated time series
c = 1 --- use the same values of t(i) for interpolation
c (t(i)=ti(i), i=1,2,...,n=nn)
c
if(is.ne.1) go to 1
if(n.ne.nn) write(*,100)
100 format(' In the polinomial fit: nn must be equal to n if is=1!')
nn=n
1 continue
c
c transform t
c
tmin=t(1)
tmax=t(n)
r=2.0d0/(tmax-tmin)
tmid=(tmin+tmax)/2.0d0
do 511 i=1,n
u(i)=t(i)
if(it.gt.0) u(i)=r*(t(i)-tmid)
511 continue
c
c change nd if the number of data is too low
c
if(n.lt.np) np=n
np1=np+1
if(np.lt.2) write(*,*) 'Too few data !'
if(np.lt.2) return
c
c compute the various power of t(i), i=1,2,...,n
c
do 500 i=1,n
td(1,i)=1.0d0
do 501 j=2,np
k=j-1
td(j,i)=u(i)*td(k,i)
501 continue
500 continue
c
c************************
c compute the fit *
c************************
c
do 1001 j=1,np
c(j)=0.0d0
1001 continue
c
c**** computation of the normal matrix
c
do 11 i=1,np
do 12 j=i,np
s=0.0d0
do 13 k=1,n
s=s+td(i,k)*td(j,k)
13 continue
g(i,j)=s
12 continue
11 continue
do 53 i=1,np
do 54 j=i,np
g(j,i)=g(i,j)
54 continue
53 continue
c
do 55 i=1,np
s=0.0d0
do 14 k=1,n
s=s+td(i,k)*x(k)
14 continue
g(i,np1)=s
55 continue
c
c**** compute the solution
c
call gelim1(g,c,np,IFLAG)
c
c compute the interpolated time series
c
dt=(tmax-tmin)/dfloat(nn-1)
c
do 520 i=1,nn
ti(i)=tmin+dt*dfloat(i-1)
if(is.eq.1) ti(i)=t(i)
time=ti(i)
if(it.gt.0) time=r*(ti(i)-tmid)
s=c(nd+1)
do 521 j=1,nd
j1=nd+1-j
s=s*time+c(j1)
521 continue
xi(i)=s
520 continue
return
end
c
c
subroutine datin8f(fname,n,t,x,bgf)
c
c-----------------------------------------------------------------------
c
c WORKS WITH F L U X E S !!
c
c-----------------------------------------------------------------------
c This routine reads in K2 time series as given by the TERRA pipeline
c by Petigure - LCs have been downloaded from the IPAC ExoFop site
c-----------------------------------------------------------------------
c
c Input: - fname = file name
c ~~~~~~
c
c Output: - {t(i),x(i),bgf(i); i=1,2,...,n}
c ~~~~~~~
c
c Comments: - Method of deriving the NORMALIZED time series:
c ~~~~~~~~~ o Employ sigma-clipping [with *avsig* clip parameter]
c to calculate the average flux
c o Divide the data by the so-obtained average flux
c - IMPORTANT: This version assumes that the comment lines
c are confined to the top *nh* lines, and ALL
c the other lines contain *NUMBERS* only, with
c proper data separators (spaces or commas)
c
c-----------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fname
character*1 hash
c
dimension t(30333),x(30333),bgf(30333)
c
c----------------------------------------------------------------------
c
nh = 40
avsig = 3.0d0
bigjd = 2400000.0d0
c
c............................................. Read in header
open(3,file=fname)
do 10 i=1,nh
read(3,*) hash
10 continue
c
c............................................. Read in data
n=0
11 continue
read(3,*,end=12) TIME_BJD,cad,xn,bgfn
n = n+1
t(n) = TIME_BJD-bigjd
x(n) = xn
bgf(n) = bgfn
go to 11
12 continue
close(3)
c
c............................................. Compute outlier-free average
call omit1(avsig,x,n,av0,si0,nclip)
c
c............................................. Normalize time series with the average
fac=1.0d0/av0
do 22 i=1,n
x(i)=x(i)*fac
22 continue
c
return
end
c
c
subroutine datin8fp(fname,n,t,x,x_p,y_p,bgf)
c
c-----------------------------------------------------------------------
c
c WORKS WITH F L U X E S !!
c
c-----------------------------------------------------------------------
c This routine reads in K2 time series as converted by "pdc2pet.f"
c from the *.tbl file resulting from the IPAC pipeline
c-----------------------------------------------------------------------
c
c Input: - fname = file name
c ~~~~~~
c
c Output: - {t(i),x(i),x_p(i),y_p(i),bgf(i); i=1,2,...,n}
c ~~~~~~~
c
c Comments: - (X,Y) pixel positions are also stored in these files
c ~~~~~~~~~ - Method of deriving the NORMALIZED time series:
c o Employ sigma-clipping [with *avsig* clip parameter]
c to calculate the average flux
c o Divide the data by the so-obtained average flux
c - IMPORTANT: This version assumes that the comment lines
c are confined to the top *nh* lines, and ALL
c the other lines contain *NUMBERS* only, with
c proper data separators (spaces or commas)
c
c-----------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fname
character*1 hash
c
dimension t(30333),x(30333)
dimension x_p(30333),y_p(30333),bgf(30333)
c
c----------------------------------------------------------------------
c
avsig = 3.0d0
bigjd = 2400000.0d0
nh = 5
c
c............................................. Read in header
open(3,file=fname)
do 10 i=1,nh
read(3,*) hash
10 continue
c
c............................................. Read in data
n=0
11 continue
read(3,*,end=12) BJD,FLUX_SAP,e_FLUX_SAP,FLUX_PDC,e_FLUX_PDC,
&X_pos,Y_pos,BG
n = n+1
t(n) = BJD-bigjd
x(n) = FLUX_SAP
x_p(n) = X_pos
y_p(n) = Y_pos
bgf(n) = BG
go to 11
12 continue
close(3)
c
c............................................. Compute outlier-free average
call omit1(avsig,x,n,av0,si0,nclip)
c
c............................................. Normalize time series with the average
fac=1.0d0/av0
do 22 i=1,n
x(i)=x(i)*fac
22 continue
c
c
return
end
c
c
subroutine read_target(itest,iseed,fname,n,t,x,iflare,
&frin0,frin1,tfreq0,dep0,tfreq1,dep1,tc0,
&x_p,y_p,ji,bgft,qtr0,qtr1,pht0,pht1,igen,std,
&xsyn,trsyn,avmag,frootv,amootv,phootv)
c
c----------------------------------------------------------------------
c Reads in TARGET time series and generates test data if itest > 0.
c----------------------------------------------------------------------
c
c Input: - itest = 0 ... Use data stored in the file
c 1 ... 2 transits + SINUSOIDAL OOTV
c 2 ... 2 transits + RR Lyrae OOTV
c 3 ... 2 transits + SINUSOIDAL OOTV + edge transits
c 4 ... 2 transits + RR Lyrae OOTV + edge transits
c 31 ... RR Lyr + 1 transit
c 32 ... RR Lyr + 1 sine
c - iseed = Big integer number to initialize random number
c generator
c - testp = Test period (in [days])
c - fname = file name
c - qtr = Fractional transit length TrTime/per,
c where 'TrTime' is the length of transit
c - dep = Depth [mag] of the transit
c - pht = Initial phase [0,1] of the transit
c - igen = 1 -- Add Gaussian noise to the noiseless
c synthetic signal
c = 0 -- Add observed time series to the
c noiseless synthetic signal
c - std = Standard deviation [mag] of the added
c Gaussian noise
c
c frin = Duration of the ingress/egress in the units
c of the total duration of the transit (from
c the first contact to the last one)
c amootv = Amplitude of the sinusoidal OOTV
c phootv = Phase of the OOTV in the units of 2*pi radian
c frootv = Frequency of the OOTV in the units of the
c orbital period
c avmag = Pre-signed average magnitude of the time series
c
c tc0 = Generated transit is assumed to have a transit
c center at: tc0 + pht/tfreq, where,(pht,tfreq) are
c are any of these: (pht0,tfreq0), (pht1,tfreq1)
c
c
c Output: - {t(i),x(i); i=1,2,...,n} = time series
c ~~~~~~~ - {tsysn(i),xsyn(i); i=1,2,...,n} = synthetic
c (noiseless) time
c series
c - ji = 1 - if we have temporal coordinates for the target
c 0 - otherwise
c
c - {x_p(i), y_p(i); i=1,2,...,n} = temporal (X,Y)
c positions for the
c target
c
c----------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fname
c
dimension t(30333),x(30333)
dimension xsyn(30333),trsyn(30333)
dimension x_p(30333),y_p(30333),bgft(30333)
c
c----------------------------------------------------------------------
c
pi2 = 8.0d0*datan(1.0d0)
testp = 1.0d0/tfreq0
c
c... Check if the target file is a derivative of an IPAC *.tbl file
c ji = 0 - it is not
c > 0 - it is
c
ji=0
ji1=index(fname,'pdc')
if(ji1.gt.0) ji=ji1
c
c... Read in data from file 'fname'
c
if(ji.eq.0) call datin8f(fname,n,t,x,bgft)
if(ji.gt.0) call datin8fp(fname,n,t,x,x_p,y_p,bgft)
do 20 i=1,n
xsyn(i)=0.0d0
trsyn(i)=0.0d0
20 continue
if(itest.eq.0) return
c
c
c***************************
c Generate test data *
c***************************
c
c........................................ The RR Lyr + 1 sine case
if(itest.eq.32) call test_32(iseed,n,t,x,iflare,
&tfreq0,dep0,tfreq1,dep1,pht2,igen,std,xsyn,avmag)
if(itest.eq.32) return
c
c........................................ The RR Lyr + 1 transit case
if(itest.eq.31) call test_31(iseed,n,t,x,iflare,
&tfreq0,tfreq1,dep1,qtr,pht,frin,igen,std,
&xsyn,trsyn,avmag)
if(itest.eq.31) return
c
c........................................ Transits with OOTV
call test_tr_oot(itest,iseed,n,t,x,iflare,
&tfreq0,tfreq1,dep0,dep1,qtr0,qtr1,pht0,pht1,tc0,
&frin0,frin1,frootv,amootv,phootv,igen,std,xsyn,trsyn)
c
return
end
c
c
subroutine test_tr_oot(itest,iseed,n,t,x,iflare,
&tfreq0,tfreq1,dep0,dep1,qtr0,qtr1,pht0,pht1,tc0,
&frin0,frin1,frootv,amootv,phootv,igen,std,xsyn,trsyn)
c
c----------------------------------------------------------------------
c This is TEST_TR_OOT (transits + OOTV + edge transits)
c----------------------------------------------------------------------
c
c NOTE: NONE
c
c======================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),trap0(30333),trap1(30333)
dimension tru(30333),xsyn(30333),trsyn(30333)
dimension x_noise(30333)
c
c----------------------------------------------------------------------
c
pi2 = 8.0d0*datan(1.0d0)
p0 = 1.0d0/tfreq0
p1 = 1.0d0/tfreq1
tim_in = 1.5d0
c................................... Parameter *tim_in* controls the epochs of
c the transits at the edges.
t14=qtr0/tfreq0
t1_1=t(1)+tim_in
t4_1=t1_1+t14
t4_2=t(n)-tim_in
t1_2=t4_2-t14
c
irr=0
if(itest.eq.2) irr=1
if(itest.eq.4) irr=1
c
c................................... Setting noise vector [with zero-average]
if(igen.eq.1) go to 3
call astat(x,n,aver,disper)
do 4 i=1,n
x_noise(i)=x(i)-aver
4 continue
go to 5
3 continue
do 6 i=1,n
x_noise(i)=std*gauss(iseed)
xsyn(i)=0.0d0
6 continue
5 continue
c
c............................................ 1st transit
e0 = tc0 + p0*pht0
call trapu(n,t,p0,e0,qtr0,frin0,tru)
do 10 i=1,n
trap0(i)=dep0*tru(i)
xsyn(i)=xsyn(i)+trap0(i)
10 continue
c
c............................................ 2nd transit
e1 = tc0 + p1*pht1
call trapu(n,t,p1,e1,qtr1,frin1,tru)
do 11 i=1,n
trap1(i)=dep1*tru(i)
xsyn(i)=xsyn(i)+trap1(i)
trsyn(i)=xsyn(i)
11 continue
c
c............................................ OOTV
do 12 i=1,n
c
e_oot = tc0 + p0*phootv
tt = pi2*(t(i)-e_oot)*frootv
xx = amootv*dsin(tt)
if(irr.eq.0) go to 99
c............................................ RR Cet (a typical RR Lyrae)
xx = 0.3193d0*dsin( tt+1.6025d0)+
&0.1558d0*dsin( 2.0d0*tt + 5.5675d0)
xx = xx + 0.1112d0*dsin(3.0d0*tt+3.6169d0)+
&0.0708d0*dsin( 4.0d0*tt + 1.6885d0)
xx = xx + 0.0408d0*dsin(5.0d0*tt+6.1338d0)+
&0.0201d0*dsin( 6.0d0*tt + 4.1837d0)
xx = xx + 0.0109d0*dsin(7.0d0*tt+2.1168d0)+
&0.0062d0*dsin( 8.0d0*tt + 0.5691d0)
xx = xx + 0.0045d0*dsin(9.0d0*tt+5.6257d0)+
&0.0038d0*dsin(10.0d0*tt + 4.4387d0)
xx=amootv*xx
c
99 continue
c
c............................................ Add transit #1 to the edges
yy=0.0d0
if(itest.eq.1) go to 15
if(itest.eq.2) go to 15
yy=0.0d0
if(t(i).lt.t1_1) go to 13
if(t(i).gt.t4_1) go to 14
yy=-dep0
go to 13
14 continue
if(t(i).gt.t4_2) go to 13
if(t(i).lt.t1_2) go to 13
yy=-dep0
13 continue
c
15 continue
xsyn(i)=xsyn(i)+xx+yy
c
12 continue
c
c............................................ Add noise to the synthetic signal
do 1 i=1,n
x(i)= xsyn(i)+x_noise(i)
1 continue
c
if(iflare.gt.0) call add_flare(n,x)
c
write(*,103) tfreq0,dep0,tfreq1,dep1
103 format(/,'This is TEST_1 (2 transits):',/,
&'f0 =',f12.7,' d0 =',f12.7,' f1 =',f12.7,' d1 =',f12.7)
c
return
end
c
c
subroutine test_31(iseed,n,t,x,iflare,
&tfreq0,tfreq1,dep1,qtr,pht,frin,igen,std,
&xsyn,trsyn,avmag)
c----------------------------------------------------------------------
c This is TEST_23 (RR Lyr + 1 transit)
c----------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension xsyn(30333),trsyn(30333)
dimension x_noise(30333)
c
c----------------------------------------------------------------------
c
call minmax(n,t,tmin,tmax)
c
pi2 = 8.0d0*datan(1.0d0)
t0 = tmin
amp1 = dep1/2.0d0
ph1 = pht
ph2 = pht+qtr
didi = frin*qtr
phin = ph1 + didi
phout = ph2 - didi
c................................... Setting noise vector [with zero-average]
if(igen.eq.1) go to 3
call astat(x,n,aver,disper)
do 4 i=1,n
x_noise(i)=x(i)-aver
4 continue
go to 5
3 continue
do 6 i=1,n
x_noise(i)=std*gauss(iseed)
6 continue
5 continue
c
do 1 i=1,n
time = t(i)-t0
c
c............................................ RR Cet
t2 = pi2*time
tt = t2*tfreq0
xx = 0.3193d0*dsin( tt+1.6025d0)+
&0.1558d0*dsin( 2.0d0*tt + 5.5675d0)
xx = xx + 0.1112d0*dsin(3.0d0*tt+3.6169d0)+
&0.0708d0*dsin( 4.0d0*tt + 1.6885d0)
xx = xx + 0.0408d0*dsin(5.0d0*tt+6.1338d0)+
&0.0201d0*dsin( 6.0d0*tt + 4.1837d0)
xx = xx + 0.0109d0*dsin(7.0d0*tt+2.1168d0)+
&0.0062d0*dsin( 8.0d0*tt + 0.5691d0)
xx = xx + 0.0045d0*dsin(9.0d0*tt+5.6257d0)+
&0.0038d0*dsin(10.0d0*tt + 4.4387d0)
c
c............................................ Transit
ph=time*tfreq1
ph=ph-idint(ph)
yy=0.0d0
rat1 = 0.0d0
rat2 = 0.0d0
if(didi.lt.1.0d-20) go to 8
rat1 = dep1/(phin-ph1)
rat2 = -dep1/(ph2-phout)
8 continue
if(ph.lt.ph1) go to 7
if(ph.gt.ph2) go to 7
yy = dep1
if(ph.lt.phin) yy = rat1*(ph-ph1)
if(ph.gt.phout) yy = dep1+rat2*(ph-phout)
7 continue
c
c............................................
zz = avmag + xx + yy
xsyn(i) = zz
trsyn(i)= yy
x(i) = zz+x_noise(i)
c
1 continue
c
if(iflare.gt.0) call add_flare(n,x)
c
write(*,103) tfreq0,tfreq1,dep1
103 format(/,'This is TEST_31 (RR Lyrae + 1 transit):',/,
&' ... Fourier decomposition of RR Cet is used ...',/,
&'tfreq0 = ',f12.7,/,
&'tfreq1 = ',f12.7,/,
&'dep1 = ',f12.7)
c
return
end
c
c
subroutine test_32(iseed,n,t,x,iflare,
&tfreq0,dep0,tfreq1,dep1,pht2,igen,std,xsyn,avmag)
c----------------------------------------------------------------------
c This is TEST_23 (RR Lyr + 1 sine)
c----------------------------------------------------------------------
c
c NOTES: - "dep0" is the factor by which the total amplitude of the
c of the prime signal is multiplied
c ... So, it is NOT the amplitude of the prime component!!
c
c - We use the light curve of RR Cet (an RR Lyrae star) with
c fundamental frequency given by the parameter "tfreq0"
c amplitude scaled by dep0.[The true period of RR_Cet is
c 0.553028770 days]
c
c - The side lobe parameters are given by (tfreq1,dep1,pht2)
c
c----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension xsyn(30333)
dimension x_noise(30333)
c
c----------------------------------------------------------------------
c
call minmax(n,t,tmin,tmax)
c
pi2 = 8.0d0*datan(1.0d0)
t0 = 0.5d0*(tmin+tmax)
amp0= dep0
amp1= dep1
ph1 = pht2
c
c................................... Setting noise vector [with zero-average]
if(igen.eq.1) go to 3
call astat(x,n,aver,disper)
do 4 i=1,n
x_noise(i)=x(i)-aver
4 continue
go to 5
3 continue
do 6 i=1,n
x_noise(i)=std*gauss(iseed)
6 continue
5 continue
c
c................................... Generate NOISELESS SYNTHETIC MONOPERIODIC time series
freq=pi2*tfreq0
do 1 i=1,n
time = t(i)-t0
tt = time*freq
xx = 0.3193d0*dsin( tt+1.6025d0)+
&0.1558d0*dsin( 2.0d0*tt + 5.5675d0)
xx = xx + 0.1112d0*dsin(3.0d0*tt+3.6169d0)+
&0.0708d0*dsin( 4.0d0*tt + 1.6885d0)
xx = xx + 0.0408d0*dsin(5.0d0*tt+6.1338d0)+
&0.0201d0*dsin( 6.0d0*tt + 4.1837d0)
xx = xx + 0.0109d0*dsin(7.0d0*tt+2.1168d0)+
&0.0062d0*dsin( 8.0d0*tt + 0.5691d0)
xx = xx + 0.0045d0*dsin(9.0d0*tt+5.6257d0)+
&0.0038d0*dsin(10.0d0*tt + 4.4387d0)
xsyn(i) = xx
1 continue
c
c................................... Find (max,min) for {xsyn}
xmin= 1.0d30
xmax=-1.0d30
pmax=0.0d0
pmin=0.0d0
do 2 i=1,n
xx=xsyn(i)
if(xx.lt.xmax) go to 7
xmax=xx
pmax=t(i)
7 continue
if(xx.gt.xmin) go to 2
xmin=xx
pmin=t(i)
2 continue
a0=xmax-xmin
t0_mod=pmax-t0+ph1/tfreq0
c
c................................... Rescale amplitude, add side lobe and noise
afac=amp0/a0
freq=pi2*tfreq1
do 8 i=1,n
time = t(i)-t0_mod
tt = time*freq
c
yy = amp1*dsin(tt)
c
xsyn(i)=afac*xsyn(i)*(1.0d0+yy)
x(i)=xsyn(i)+x_noise(i)+avmag
8 continue
c
if(iflare.gt.0) call add_flare(n,x)
c
write(*,103) tfreq0,tfreq1,amp1
103 format(/,'This is TEST_32 (RR Lyrae + 1 sine):',/,
&' ... Fourier decomposition of RR Cet is used ...',/,
&'tfreq0 = ',f12.7,/,
&'tfreq1 = ',f12.7,/,
&'amp1 = ',f12.7)
c
return
end
c
c
subroutine read_ref(n_ref,ref_id,ref_ra,ref_de,ref_mag)
c---------------------------------------------------------------------
c This routine reads in reference file
c---------------------------------------------------------------------
c
c INPUT - ref.pos = Reference/position/magnitude file
c
c OUTPUT: - n_ref = Number of reference stars
c - {ref_id} = IDs of the reference stars
c - {ref_ra} = RAs of the reference stars
c - {ref_de} = DEs of the reference stars
c - {ref_mag} = Average magnitudes of the reference stars
c
c NOTES: - NONE
c
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 ref_id(30333)
c
dimension ref_mag(30333),ref_ra(30333),ref_de(30333)
c
c--------------------------------------------------------------------
c
c... Get the number of the available targets
c
call SYSTEM('wc -l ref.pos > nr.dat')
open(1,file='nr.dat')
read(1,*) n_ref
close(1)
c
c... Read in reference magnitudes
c
open(1,file='ref.pos')
do 5 i=1,n_ref
read(1,*) ref_id(i),ref_ra(i),ref_de(i),ref_mag(i)
5 continue
close(1)
c
return
end
c
c
subroutine file_in(fname,sname,nfile,camp)
c---------------------------------------------------------------------------
c This routine reads in the next star name, path and K2 campaing number
c---------------------------------------------------------------------------
c
c INPUT: - filename.lis
c ** Per line format: 211681517 LC/211681517-c05.txt
c
c OUTPUT: - fname : Full path to the time series
c - sname : Star name (9 characters)
c - nfile : Number of the files remained in 'filename.lis'
c - camp : Campampaing ID (e.g., c03)
c
c NOTES: - After each call of this routine the first row in 'filename.lis'
c will be deleted ...
c - The file format is *strict* as given above, i.e.:
c * star name is a 9 character string
c * campaign number must be in the 14-16 character range of
c the full path
c * the full path cannot be greater than 25
c
c---------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*1 blank
character*3 camp
character*9 sname,sn(30333)
character*25 fname,path(30333)
c
c---------------------------------------------------------------------------
blank =' '
c
open(13,file='filename.lis')
c
read(13,7,end=15) sname,fname
7 format(a9,9x,a25)
15 continue
if(sname.eq.blank) nfile=-1
if(nfile.eq.-1) write(*,100)
100 format('No more file names in *filename.lis* !')
if(sname.eq.blank) return
camp = fname(14:16)
c
c... Save the remaining content of the file
c
n = 0
8 continue
n = n+1
read(13,7,end=3) sn(n),path(n)
go to 8
3 continue
n = n-1
close(13)
c
open(13,file='filename.lis')
if(n.eq.0) go to 103
do 10 i=1,n
write(13,7) sn(i),path(i)
10 continue
go to 102
103 continue
write(13,7) blank
102 continue
close(13)
nfile=n
c
return
end
c
c
subroutine sppolfit_w(ndsp,nf,p)
c
c-------------------------------------------------------------------
c This routine fits a polynomial to the frequency spectrum
c and returns the residual spectrum. Weighted LS is employed
c-------------------------------------------------------------------
c
c INPUT: - ndsp = degree of the polynomial to be fitted to the
c spectrum
c - nf = number of the spectrum points
c - p(i) = spectrum array
c
c OUTPUT: p(i) = spectrum array [trend-free]
c
c NOTES: - It is assumed that p(i) is equidistant with
c respect of the frequency.
c - For saving time, an evenly sampled subset of the
c full spectrum is fitted.
c
c-------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),p(155555)
dimension c(995)
c
if(ndsp.lt.1) return
c
nd = ndsp
nsm = 30333
c
c nsm = Upper limit of the number of spectrum points used in
c the polynomial fit. These points are distributed
c uniformly in the frequency band.
c-------------------------------------------------------------------
c
n2 = 1+idint(dfloat(nf)/dfloat(nsm))
n = 0
do 4 i=1,nf,n2
n = n+1
t(n) = dfloat(i)
x(n) = p(i)
4 continue
c
r=2.0d0/(t(n)-t(1))
tmid=(t(1)+t(n))/2.0d0
c
call robust_pol2(nd,n,t,x,c,sig)
c
c Subtract fitted polynomial
c
pmin = 1.0d30
do 1 i=1,nf
s=c(nd+1)
time=dfloat(i)
time=r*(time-tmid)
do 2 j=1,nd
j1=nd+1-j
s=s*time+c(j1)
2 continue
p(i) = p(i)-s
pmin = dmin1(pmin,p(i))
c
1 continue
c
do 3 i=1,nf
p(i) = p(i)-pmin
3 continue
c
return
end
c
c
subroutine read_tfa_lis(ji,nte,fid,idt0)
c-----------------------------------------------------------
c This routine reads in pre-defined list of templates
c-----------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fname,fid(995)
character*20 fn
character*9 sname,idt0(995)
character*2 lc
character*1 slash
c
c... Check data type on *template.lis*
c
k=0
open(30,file='template.lis')
read(30,2,end=3) sname,lc,slash,fn
k=k+1
2 format(a9,9x,a2,a1,a20)
fname=trim(lc//slash//fn)
3 continue
close(30)
if(k.eq.0) write(*,4)
if(k.eq.0) stop
4 format('>>>>>> template.lis is EMPTY !')
c
c... Check if TFA template star names are consistent with the
c data type of the target
c
k=index(fname,'pdc')
if(k.ne.ji) write(*,1)
if(k.ne.ji) stop
1 format(56('#'),/,
&' TFA template data type contradicts to TARGET data type',/,
&56('#'))
c
if(k.gt.0) go to 5
c
c....................................... Read in TERRA-type list
nte=0
open(30,file='template.lis')
600 continue
nte=nte+1
read(30,603,end=601) idt0(nte),fid(nte),ih
603 format(a9,9x,a20,1x,i4)
go to 600
601 continue
nte=nte-1
close(30)
go to 800
c
c....................................... Read in PDC-type list
5 continue
nte=0
open(30,file='template.lis')
700 continue
nte=nte+1
read(30,703,end=701) idt0(nte),fid(nte),ih
703 format(a9,9x,a23,1x,i5)
go to 700
701 continue
nte=nte-1
close(30)
c
800 continue
c
return
end
c
c
subroutine tfa_ts(itemp,ji,nte,fid,mte0,idt0)
c---------------------------------------------------------------------
c This routine reads in TFA, BG and XY time series
c---------------------------------------------------------------------
c
c INPUT: - itemp : 0 - if no TFA template is to be used
c >0 - if the TFA templates given in the
c template file list is to be used
c - ji : 0 - ExoFOP download (Petigura's TERRA)
c : 1 - IPAC download (Kepler PDC pipeline)
c - nte : Number of templates
c (i.e., number of array elements of {fid})
c - {fid} : Input file IDs
c - {idt0} : Input object names
c
c OUTPUT: - {xte0} : Template time series
c -->common /aa0/
c - {tte0} : Moments of time of {xte0}
c --> common /aa0/
c - {nte0} : Number of data points of {xte0}
c --> common /aa0/
c - {idt0} : Output object names [see NOTES]
c - mte0 : nte+1
c
c NOTE: - After completion: nte --> nte+1
c [the constant template is added]
c - After completion: The object name array {idt0} is
c shifted by one to allow constant
c template ID.
c
c=====================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fid(995)
character*25 fname
character*9 idt(995),idt0(995)
c
dimension t(30333),x(30333),bgf(30333),x_p(30333),y_p(30333)
c
common /aa0/ tte0(995,30333),xte0(995,30333),nte0(995)
common /tt/ ttemp(30333),ntemp
c
c................................... Set CONSTANT template and the current
c template timebase
k=1
do 599 i=1,ntemp
xte0(k,i) = 1.0d0
599 continue
nte0(k)=ntemp
idt(k)='111111111'
if(itemp.eq.0) go to 200
c
do 4 kk=1,nte
c
c................................... Read in template
fname = fid(kk)
if(ji.eq.0) call datin8f(fname,n,t,x,bgf)
if(ji.gt.0) call datin8fp(fname,n,t,x,x_p,y_p,bgf)
k=k+1
idt(k)=idt0(kk)
nte0(k)=n
c
c................................... Store the corresponding template array
do 44 i=1,n
tte0(k,i)=t(i)
xte0(k,i)=x(i)
44 continue
c
4 continue
c
200 continue
c................................... Store template names
do 201 j=1,k
idt0(j)=idt(j)
201 continue
mte0=k
c
return
end
c
c
subroutine temp_new(nseq,ji,idvar,itemp,x_p,y_p,bgft,ibg,ixy,
&nt,ta,xa,intemp,x_int)
c
c**************************************************************************
c >>> Read in / reload TFA template, BG and (X,Y) time series
c >>> Interpolate the above time series to the timebase of the target
c**************************************************************************
c
c INPUT: nseq - 1 for the FIRST call of this routine
c >1 for ALL SUBSEQUENT calls
c ji - 0 - ExoFOP download is used (based on TERRA)
c 1 - IPAC download is used (based on the Kepler pipe)
c idvar - Character ID of the variable to be analyzed
c itemp - 1 use TFA templates for co-trending
c 0 do NOT use TFA templates
c {x_p} - Temporal X coordinate of the target
c {y_p} - Temporal Y coordinate of the target
c {bgft} - Bacground flux of the target
c ibg - Key for target background flux for LC correction
c ixy - Key for target position for LC correction
c nt - Number of data points of the target time series
c ta(i) - Time values of the target time series
c xa(i) - Values of the target time series
c
c
c OUTPUT: x_int(i) - Input time series (i.e., {x_int}={xa})
c xte(j,i) - Current template set ( common /aa/ )
c xte0(j,i)- Template set of the 1st call ( common /aa0/ )
c tte0(j,i)- Original timebase of {xte0} ( common /aa0/ )
c mte0 - TFA template number (without BG and XY vectors)
c set in the 1st call
c m3_tfa - TFA+(BG, XY vector) [because of the avaibility
c of the BG, XY vectors, m3_tfa is specific to
c each call]
c
c
c NOTES: - NO CHANGE in the input time series {ta,xa}
c - TFA template timebase = Target timebase
c - The array {x_int} comes from an earlier version of this
c routine (when we used a common timebase not equivalent with
c that of the target - so, we had to interpolate the target
c time series to the template timebase - this is why we have
c "int" in x_int). This has been left in, to avoid changes in
c other places.
c
c--------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*25 fid(995)
character*9 idvar
character*9 idt(995),idt0(995)
c
dimension t(30333),x(30333),y(30333),ta(30333),xa(30333)
dimension x_int(30333)
dimension x_p(30333),y_p(30333),bgft(30333)
c
common /aa/ xte(995,30333)
common /aa0/ tte0(995,30333),xte0(995,30333),nte0(995)
common /cc/ m3_tfa,mte0
common /cc1/ idt
common /tt/ ttemp(30333),ntemp
c
c--------------------------------------------------------------------------
c
c SET BASIC PARAMETERS
c
nbg = 3
nxy = 3
intemp= 0
c
c nbg = Degree of the polynomial of the background detrending
c nxy = Degree of the polynomial of the position detrending
c
c............................................ Set the template timebase
c equal to the target timebase
c and {x_int}={xa}
c
ntemp=nt
do 4 i=1,ntemp
ttemp(i)=ta(i)
x_int(i)=xa(i)
xte(1,i)=1.0d0
4 continue
c
c.... Read in pre-defined list of templates
c
if(itemp.eq.1.and.nseq.eq.1) call read_tfa_lis(ji,nte,fid,idt0)
c
c.... Read in template time series
c
if(nseq.eq.1) call tfa_ts(itemp,ji,nte,fid,mte0,idt0)
mte1=mte0
if(itemp.eq.0) go to 6
c
c... We have the basic template set in the initialization part above:
c
c {nte0(j),tte0(j,i),xte0(j,i); j=1,2,...,mte0; i=1,2,...,nte0(j)}
c
c Here we: 1. INTERPOLATE TFA template set to the timebase of the
c target
c 2. CHECK if the target is a member of the template set
c 3. ADD other systematics correction vectors
c
c.................................... Interpolate TFA template to {ttemp}
c
do 11 j=2,mte0
n=nte0(j)
do 13 i=1,n
t(i)=tte0(j,i)
x(i)=xte0(j,i)
13 continue
call tbase_inter(n,t,x,y)
do 12 i=1,ntemp
xte(j,i)= y(i)
12 continue
11 continue
c.................................... Check if the target is a member
c of the template set. If yes, then
c exclude the corresponding template
it=0
do 5 i=2,mte0
if(idt0(i).eq.idvar) it=i
5 continue
if(it.eq.0) go to 6
c
c.................................... Exclude the template with the same
c identity as the target
intemp=1
mte1=mte0-1
if(it.eq.mte0) go to 6
do 7 j=it,mte1
j1=j+1
do 8 i=1,ntemp
xte(j,i)=xte(j1,i)
8 continue
7 continue
c
6 continue
c
c............................................. Add BACGROUND vector to the templates
c [a nbg-order polynomial]
kin=mte1
k=kin
if(ibg.gt.0) call add_bgf(kin,nbg,bgft,kout)
if(ibg.gt.0) k=kout
c
c............................................. Add POSITION vectors to the templates
c [a nxy-order polynomial]
c
kin=k
if(ixy.gt.0) call add_xy(kin,nxy,x_p,y_p,kout)
if(ixy.gt.0) k=kout
c
nte=k
write(*,144) nte
144 format('Total number of TFA/EPD correcting vectors =',i4)
c
m3_tfa=nte
c
return
end
c
c
subroutine add_bgf(kin,nbg,bgft,kout)
c-----------------------------------------------------------
c Add BACGROUND vector to the templates
c [powers of the background flux]
c-----------------------------------------------------------
c
c INPUT:
c kin = Number of TFA templates (including constant)
c at the moment of the call of this routine
c nbg = Degree of the polynomial of the background
c detrending vector (each polynomial order
c increases the number of correcting vectors
c associated with the background)
c {bgft}= Background flux time series
c
c OUTPUT:
c kout = Total number of co-trending vectors at the
c moment of exiting from this routine
c {xte}= Updated TFA vectors with the background
c vectors
c
c===========================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension bgft(30333),b(30333),t(30333),y(30333)
c
common /aa/ xte(995,30333)
common /tt/ ttemp(30333),ntemp
c
c-----------------------------------------------------------
c
k=kin
eps=1.0d-10
np=20
asig3=3.0d0
c
c eps = if the RMS of the background flux is greater than
c eps, then NO background flux is added to the TFA set
c np = order of the polynomial to be fitted to the
c background flux
c asig3= sigma clipping parameter for the background flux
c
c... Compute relative background flux vector
c
call astat(bgft,ntemp,av_bg,sigma)
if(sigma.gt.eps) go to 4
nbg=0
go to 6
4 continue
do 10 i=1,ntemp
t(i) = ttemp(i)
b(i) = bgft(i)
10 continue
c
c... Treating outliers based on a *np*-order polynomial fit
c
ipr=-1
n=ntemp
call robust_pol(n,np,t,b,y,sig,stdv)
call sigma_clipping_r(ipr,asig3,n,b,y,nclip)
c
c... Compute powers of {b}
c
do 1 m=1,nbg
k=k+1
do 3 i=1,ntemp
xx=b(i)
sx=1.0d0
do 5 j=1,m
sx=sx*xx
5 continue
y(i)=sx
3 continue
call astat(y,ntemp,aver,sigma)
fac=1.0d0/aver
do 2 i=1,ntemp
xte(k,i)=fac*y(i)-1.0d0
2 continue
1 continue
c
6 continue
kout=k
write(*,558) kin,nbg,kout
558 format(/,'>>> Adding BACKGROUND vectors:',/,
& 'Starting number of templates+1 =',i4,/,
& 'Number of background vectors =',i4,/,
& 'Final number of templates =',i4,/)
c
return
end
c
c
subroutine add_xy(kin,nxy,x_p,y_p,kout)
c-----------------------------------------------------------
c Add POSITION vectors to the templates
c [powers of the (X,Y) positions]
c-----------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),y(30333),z(30333)
dimension x_p(30333),y_p(30333)
c
common /aa/ xte(995,30333)
common /tt/ ttemp(30333),ntemp
c
c-----------------------------------------------------------
c
k=kin
kout=kin
if(nxy.eq.0) return
c
c... Compute relative (X,Y) vectors
c
call astat(x_p,ntemp,avx_p,sigma)
call astat(y_p,ntemp,avy_p,sigma)
avx_p=1.0d0/avx_p
avy_p=1.0d0/avy_p
do 260 i=1,ntemp
x(i) = x_p(i)*avx_p-1.0d0
y(i) = y_p(i)*avy_p-1.0d0
260 continue
c
c... Compute powers of {x},{y}
c
kxy=0
do 1 m=1,nxy
m1=m+1
do 2 j=1,m1
kxy=kxy+1
k=k+1
jx=j-1
jy=m1-j
do 3 i=1,ntemp
xx=x(i)
yy=y(i)
c.................................. X**jx
sx=1.0d0
if(jx.eq.0) go to 6
do 5 jp=1,jx
sx=sx*xx
5 continue
6 continue
c.................................. Y**jy
sy=1.0d0
if(jy.eq.0) go to 8
do 7 jp=1,jy
sy=sy*yy
7 continue
8 continue
z(i)=sx*sy
3 continue
call astat(z,ntemp,aver,sigma)
fac=1.0d0/aver
do 9 i=1,ntemp
xte(k,i)=fac*z(i)-1.0d0
9 continue
2 continue
1 continue
c
kout=k
write(*,558) kin,kxy,kout
558 format('>>> Adding POSITION vectors:',/,
& 'Starting number of templates+1 =',i4,/,
& 'Number of position vectors =',i4,/,
& 'Final number of templates =',i4,/)
c
return
end
c
c
subroutine tbase_inter(n,t,x,y)
c-----------------------------------------------------------------
c This routine linearly interpolates time series
c {x(i); i=1,2,...,n} to the timebase given by
c {ttemp(i); i=1,2,...,ntemp}
c Output time series: {y(i); i=1,2,...,n}
c
c NOTES: - Returns array {y(i); i=1,2,...,ntemp} that is:
c * zero-averaged
c * defined on the template timebase
c {ttemp(i); i=1,2,...,ntemp}
c
c - It is assumed that both time bases - {t} and {temp}
c are monotinic
c-----------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 idt(995)
c
common /cc/ m3_tfa,mte0
common /cc1/ idt
common /tt/ ttemp(30333),ntemp
c
dimension t(30333),x(30333),y(30333)
c
c.............................. FIRST-order interpolation
c NOTE: If extrapolation is needed, we set
c the values equal to the first/last
c value of the time series.
c
c... Compute average of {x}
c
s=0.0d0
do 6 i=1,n
s=s+x(i)
6 continue
aver=s/dfloat(n)
c
s=0.0d0
xn=x(n)
do 1 i=1,ntemp
tt = ttemp(i)
k = 0
y(i) = aver
do 2 j=1,n
if(k.gt.0) go to 2
if(t(j).lt.tt) go to 2
k=j
2 continue
if(k.eq.0) go to 3
if(k.eq.1) k=2
k2=k
k1=k2-1
x1=t(k1)
x2=t(k2)
y1=x(k1)
y2=x(k2)
yy = y1 + (tt-x1)*(y2-y1)/(x2-x1)
go to 5
3 continue
yy=xn
5 continue
s=s+yy
y(i)=yy
1 continue
s=s/dfloat(ntemp)
c
do 4 i=1,ntemp
y(i)=y(i)-s
4 continue
c
return
end
c
c
subroutine spsnr1(nf,p,fmin,df,frmax,pmax,snr,spd)
c---------------------------------------------------------------------
c This routine computes the Signal-to-Noise Ratio and the Spectral
c Peak Density of the frequency spectrum
c---------------------------------------------------------------------
c
c INPUT: nf = number of points in the frequency spectrum
c p(i) = frequency spectrum , i=1,2,...,nf
c fmin = frequency at p(1)
c df = frequency step
c
c OUTPUT: frmax = Frequency at the maximum peak
c pmax = Spectrum value at the maximum peak
c snr = Signal-to-Noise Ratio
c spd = Spectral Peak Density
c
c NOTE: - Signal-to-Noise Ratio is defined as follows:
c [ p(highest peak) - ] / sigma(p)
c
c - If asig6 > 666.0, then sigma clipping is skipped for
c the calculation of
and sigma(p)
c
c - SPD was introduced on 2020.04.25 to make distinction
c between SPARSE and DENSE spectra. This parameter is
c equal to the ratio of the frequency regime occupied
c by spectrum values above *r_peak*, the fraction of
c the amplitude of the highest peak, above which the
c occupation rate is calculated.
c
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension p(155555),y(155555)
dimension in(155555)
c
asig6 = 777.0d0
r_peak= 0.30d0
eps0 = 1.0d-15
nn = nf
c
c................................. Find maximum of the spectrum
pmax=-1.0d30
k=0
do 10 i=1,nn
if(p(i).lt.pmax) go to 10
pmax=p(i)
k=i
10 continue
frmax=fmin+df*dfloat(k-1)
c
if(asig6.lt.666.0d0) go to 13
c
c................................. Compute RMS of {p} without sigma clipping
s=0.0d0
do 14 i=1,nn
s=s+p(i)
14 continue
s=s/dfloat(nn)
d=0.0d0
do 15 i=1,nn
dd=p(i)-s
d=d+dd*dd
15 continue
xave=s
xsig=dsqrt(d/dfloat(nn-1))
go to 40
c
13 continue
c
c................................. Store working array
do 12 i=1,nn
y(i)=p(i)
12 continue
20 continue
c
c................................. Compute average, sigma
s=0.0d0
do 30 i=1,nn
s=s+y(i)
30 continue
xave=s/dfloat(nn)
dev=0.0d0
do 31 i=1,nn
d=y(i)-xave
dev=dev+d*d
31 continue
xsig=dsqrt(dev/dfloat(nn-1))
if(xsig.le.eps0) go to 40
c
c................................. Find outliers
devm=asig6*xsig
k=0
do 11 i=1,nn
dev=dabs(y(i)-xave)
j=1
if(dev.gt.devm) j=0
k=k+j
in(i)=j
11 continue
c................................. If no outlier found, compute SDE
if(k.eq.nn) go to 22
c
c................................. Omit outliers
k=0
do 1 i=1,nn
if(in(i).eq.0) go to 1
k=k+1
y(k)=y(i)
1 continue
nn=k
go to 20
c
22 continue
c
c................................. Sigma clipping completed,
c update average and sigma
s=0.0d0
do 32 i=1,nn
s=s+y(i)
32 continue
xave=s/dfloat(nn)
dev=0.0d0
do 33 i=1,nn
d=y(i)-xave
dev=dev+d*d
33 continue
xsig=dsqrt(dev/dfloat(nn-1))
c
40 continue
c
snr = (pmax-xave)/xsig
if(snr.gt.99.9d0) snr=99.9d0
c
c................................. Compute SPD
c
np=0
rpmax=1.0d0/pmax
do 50 i=1,nf
r=p(i)*rpmax
if(r.gt.r_peak) np=np+1
50 continue
spd=dfloat(np)/dfloat(nf)
c
return
end
c
c
subroutine omit1(asig,x,n,xave,xsig,nclip)
c
c---------------------------------------------------------------------
c
c This routine computes average and standard deviation
c of { x(i) ; i = 1,2, ..., n } by iteratively omitting
c data which deviate more that asig*sigma
c
c Input parameters/arrays: asig, t(i), x(i), n
c ~~~~~~~~~~~~~~~~~~~~~~~~
c
c Output parameters: xave = average
c ~~~~~~~~~~~~~~~~~~ xsig = standard deviation
c nclip = number of data points clipped
c
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),y(30333)
dimension in(30333)
c
eps0 = 1.0d-15
nn = n
c
do 21 i=1,nn
y(i) = x(i)
21 continue
c
20 continue
c
call astat(y,nn,aver,disper)
c
c... Finding out if there is an outlier
c
if(disper.le.eps0) return
c
iout=0
devm=asig*disper
do 10 i=1,nn
dev=dabs(y(i)-aver)
in(i)=1
if(dev.gt.devm) in(i)=0
10 continue
c
c... Omit outlier
c
k=0
do 1 i=1,nn
if(in(i).eq.0) go to 1
k=k+1
y(k)=y(i)
1 continue
if(k.eq.nn) go to 22
nn=k
go to 20
22 continue
nclip=n-nn
call astat(y,nn,xave,xsig)
c
return
end
c
c
subroutine tru_fit(n,trap,x,c0,c1,sdev)
c
c-------------------------------------------------------------------
c This routine performs the following LS fit:
c
c SUM [x3(i) - c0 - c1*trap(i)]**2 = minimum
c
c Input: n = Number of data points
c ~~~~~~ trap(i) = Function to be fitted to the target function
c x(i) = Target funcion
c
c Output: c0, c1 = Regression coefficients
c ~~~~~~~ sdev = Standard deviation of the residuals of
c the fit [biased estimate]
c
c Notes: - The two time series MUST be defined on the SAME
c time base, no empty values in both of them
c
c-------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension trap(30333),x(30333)
c
n0=n
maxdim = 30333
if(n.gt.maxdim) write(*,*) ' n > maxdim in *tru_fit* !'
if(n.gt.maxdim) n=maxdim
c
c... Compute normal matrix
c
a11=0.0d0
a12=0.0d0
a22=0.0d0
a13=0.0d0
a23=0.0d0
do 1 i=1,n
a=trap(i)
b=x(i)
a12=a12+a
a22=a22+a*a
a13=a13+b
a23=a23+b*a
1 continue
a11=dfloat(n)
a21=a12
d=1.0d0/(a11*a22-a12*a21)
c0 =(a13*a22-a23*a12)*d
c1 =(a23*a11-a13*a21)*d
c
s=0.0d0
do 2 i=1,n
d=x(i)-c0-c1*trap(i)
s=s+d*d
2 continue
sdev = dsqrt(s/dfloat(n))
n=n0
c
return
end
c
c
subroutine tru_fit_w(n,we,trap,x,c0,c1,sdev)
c
c-------------------------------------------------------------------
c This routine performs the following WEIGHTED LS fit:
c
c SUM we(i)*[x3(i) - c0 - c1*trap(i)]**2 = minimum
c
c Input: n = Number of data points
c ~~~~~~ we(i) = Point-by-point weights
c trap(i) = Function to be fitted to the target function
c x(i) = Target funcion
c
c Output: c0, c1 = Regression coefficients
c ~~~~~~~ sdev = Standard deviation of the residuals of
c the fit [biased estimate]
c
c Notes: - The two time series MUST be defined on the SAME
c time base, no empty values in both of them
c
c-------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension trap(30333),x(30333),we(30333)
c
n0=n
maxdim = 30333
if(n.gt.maxdim) write(*,*) ' n > maxdim in *tru_fit* !'
if(n.gt.maxdim) n=maxdim
c
c... Compute normal matrix
c
a11=0.0d0
a12=0.0d0
a22=0.0d0
a13=0.0d0
a23=0.0d0
do 1 i=1,n
w=we(i)
a=trap(i)
b=x(i)
a11=a11+w
a12=a12+w*a
a22=a22+w*a*a
a13=a13+w*b
a23=a23+w*b*a
1 continue
a21=a12
d=1.0d0/(a11*a22-a12*a21)
c0 =(a13*a22-a23*a12)*d
c1 =(a23*a11-a13*a21)*d
c
s=0.0d0
do 2 i=1,n
d=x(i)-c0-c1*trap(i)
s=s+d*d
2 continue
sdev = dsqrt(s/dfloat(n))
n=n0
c
return
end
c
c
subroutine events(n,t,f0,epoch,qtran,nev,ndit,itoc)
c
c========================================================================
c This routine computes number of transits with observed
c data points in them
c
c Input: n -- Number of data points
c ~~~~~~ t(i) -- Time values
c f0 -- Folding frequency (1/period)
c epoch -- Epoch of the first ingress
c qtran -- Fractional transit length
c (in the units of the period)
c
c Output: nev -- Number of transits with observed data points
c ~~~~~~~ in them
c ndit -- Total number of data points in the transit
c itoc -- Transit Occupation Code
c An *ntc* digit integer, giving the relative
c number of data points in the top occupied
c transit events. The occupation ratio decreases
c from left to right. The highest occupation
c rate, corresponding to 10, is not shown. The
c leftmost digit shows the ratio of the number
c of data points of the next most occupated
c transit event to the most occupated one.
c That is, if this digit is, e.g., 7, this means
c that the second most occupied event contains
c 0.7*N_max data points, where N_max is the
c number of data points sitting in the event
c containing the highest number of data points.
c========================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),tr(30333)
c
c========================================================================
c
nevm = 30333
ntc = 5
p0 = 1.0d0/f0
dt = p0*qtran
call minmax(n,t,tmin,tmax)
c
c... Find the first moment of transit preceeding *tmin*
c
t0=epoch
11 continue
if(t0.lt.tmin) go to 12
t0=t0-p0
go to 11
12 continue
c
c... Find number of transits with observed data points
c
ndit=0
nev=0
ne=-1
trm=-1.0d30
k=0
1 continue
ne=ne+1
t1=t0+p0*dfloat(ne)
if(t1.gt.tmax) go to 10
t2=t1+dt
k=0
do 2 i=1,n
time=t(i)
if(time.lt.t1) go to 2
if(time.gt.t2) go to 2
k=k+1
2 continue
if(k.lt.1) go to 1
if(nev.gt.nevm) go to 10
nev=nev+1
ndit=ndit+k
tr(nev)=dfloat(k)
trm=dmax1(trm,tr(nev))
go to 1
10 continue
c
c... Compute Transit Occupation Code [TOC]
c
ntc1=ntc+1
toc=0.0
do 3 j=1,ntc1
tcm=-1.0d30
do 4 i=1,nev
if(tr(i).lt.tcm) go to 4
k=i
tcm=tr(i)
4 continue
tr(k)=-1.0
if(j.eq.1) go to 3
if(tcm.lt.0.0) tcm=0.0
tcm=idint(10.0*tcm/trm)
tcm=dmin1(9.0d0,tcm)
toc=toc+tcm*(10.0**(ntc1-j))
3 continue
itoc=idint(toc)
c
return
end
c
C
SUBROUTINE INTER2(IQ,X,Y,N1,XD,YD,N2)
C
C-----------------------------------------------------------------
C THIS ROUTINE INTERPOLATES THE SEQUENCE
C X(I),Y(I), I=1,2,...,N1
C
C IN X(I), INTO THE OTHER GRID GIVEN BY {XD(I), I=1,2,...,N2}
C
C OUTPUT IS GIVEN IN:
C XD(I),YD(I), I=1,2,...,N2
C
C VERSION: 2010.04.01 (cleaned up from INTERX)
C-----------------------------------------------------------------
C
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
C
DIMENSION Y(30333),X(30333)
DIMENSION XD(30333),YD(30333)
C
C IQ = 2 ... QUADRATIC INTERPOLATION
C = 1 ... LINEAR INTERPOLATION
C
DO 10 I=1,N2
C
C FIND THE PROPER ARRAY INDICES FOR THE INTERPOLATION REGIME
C
RX=XD(I)
DEV=1.0D30
JJ=1
DO 1 J=1,N1
DD=DABS(RX-X(J))
IF(DD.GT.DEV) GO TO 1
DEV=DD
JJ=J
1 CONTINUE
IF(JJ.EQ.1) GO TO 2
IF(JJ.EQ.N1) GO TO 3
J1=JJ-1
J2=JJ
J3=JJ+1
GO TO 4
2 CONTINUE
J1=1
J2=2
J3=3
GO TO 4
3 CONTINUE
J1=N1-2
J2=N1-1
J3=N1
4 CONTINUE
C
IF(IQ.EQ.1) GO TO 5
C
C QUADRATIC INTERPOLATION
C
AA=(RX-X(J2))*(RX-X(J3))/((X(J1)-X(J2))*(X(J1)-X(J3)))
BB=(RX-X(J1))*(RX-X(J3))/((X(J2)-X(J1))*(X(J2)-X(J3)))
CC=(RX-X(J1))*(RX-X(J2))/((X(J3)-X(J1))*(X(J3)-X(J2)))
YD(I)=Y(J1)*AA+Y(J2)*BB+Y(J3)*CC
GO TO 11
C
5 CONTINUE
C
C LINEAR INTERPOLATION
C
J1=JJ
IF(J1.EQ.1) GO TO 6
IF(J1.EQ.N1) GO TO 7
J2=J1+1
IF(DABS(RX-X(J1-1)).LT.DABS(RX-X(J1+1))) J2=J1-1
GO TO 8
6 CONTINUE
J2=2
GO TO 8
7 CONTINUE
J2=N1-1
8 CONTINUE
C
YD(I)=Y(J1)+(RX-X(J1))*(Y(J2)-Y(J1))/(X(J2)-X(J1))
C
11 CONTINUE
C
10 continue
C
RETURN
END
C
c
subroutine dsp_ok(n,t,x,p0,tcen,qtran,dsp,snrt,fret)
c
c===============================================================================
c This routine computes DSP and OOTV parameters SNRT and FRET
c===============================================================================
c
c INPUT:
c
c n = Number of data points
c {t} = Time values of the time series
c {x} = Values of the time series
c p0 = Period
c tcen = Moment of time of the CENTER of the TRANSIT
c qtran = Relative length of the transit (t14/P)
c
c
c OUTPUT:
c
c dsp = depth / SQRT( var_depth + var_oot + dif_oot)
c
c snrt = SNR of the highest peak in the Fourier frequency spectrum
c of the OOTV
c
c fret = Frequency [in the units of p0] at the highest peak in the
c Fourier frequency spectrum of the OOTV
c
c NOTES: - The goal in the computation of the new version of DSP is to
c consider ALL parts of the folded time series EQUALLY in
c the statistical sense. Therefore, we aim at the SAME BINNING
c as given by the TRANSIT LENGTH. The time series is phased
c in [-0.5,0.5] with the transit center at phase zero.
c End phases may not be properly covered, because of the general
c non-commensurate nature of the transit length and the full
c phase.
c
c - The scatter of the folded LC is characterized by two parameters:
c var_depth = VARIANCE of the AVERAGE of the INTRANSIT datapoints
c var_oot = VARIANCE of the OOT bin AVERAGES
c diff_oot = AVERAGE of the SQUARED DIFFERENCES of the bin
c averages in the OOT regime
c
c - The intransit variance is computed with the assumption of
c flat-bottomed shape. Therefore, deviation from this shape
c (i.e., long ingress/egress phases) lowers the final value
c of DSP [this is helpful in putting less weights on likely
c EB false positives or on other variables with V-shaped dips].
c
c - Explanation of the usage of var_oot and diff_oot:
c var_oot : whenever there is an OOTV, this parameter is in
c effect (independently of the smoothness of the
c OOT function)
c diff_oot: for ragged/discontinuous OOT, the differences
c are more relevant, than the variance. Admittedly,
c this parameter will also decrease DSP for multis,
c not only for singles with untreated outliers.
c
c - The bin variance is computed with the following formula:
c
c SI = (1/(n-1))*SUM(x_i**2) - (n /(n-1)*(**2)
c = (1/n)*SUM(x_i)
c
c===============================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),u(30333),v(30333),y(30333)
dimension p(155555)
dimension ph1(30333),ph2(30333)
dimension av(30333),bv(30333)
dimension in(30333)
c
c------------------------------------------------------------------------
c
f0=1.0d0/p0
c
c... Set bin borders with transit bin index *jc*
c
n1=idint(0.5d0/qtran-0.5d0)
nbq=2*n1+1
qtran2=qtran*0.5d0
do 29 j=1,nbq
ph1(j)=dfloat(j-n1-1)*qtran-qtran2
ph2(j)=ph1(j)+qtran
29 continue
jc=n1+1
c
c... Compute bin statistics according to the above bin distribution
c
do 28 j=1,nbq
in(j)=0
av(j)=0.0d0
28 continue
c
n3=0
njc=0
s=0.0d0
ping=ph1(jc)
pegr=ph2(jc)
c
do 30 i=1,n
ph = f0*(t(i)-tcen)
ph = ph - idint(ph)
if(ph.gt. 0.5d0) ph=ph-1.0d0
if(ph.lt.-0.5d0) ph=ph+1.0d0
xx=x(i)
if(ph.le.ping) go to 3
if(ph.ge.pegr) go to 3
njc=njc+1
y(njc)=xx
go to 4
c
c.................................... Collect OOT values
3 continue
n3=n3+1
u(n3)=ph
v(n3)=xx
s=s+xx
c
c.................................... Update bin 1st moments
4 continue
do 31 j=1,nbq
if(ph.lt.ph1(j)) go to 31
if(ph.gt.ph2(j)) go to 31
in(j)=in(j)+1
av(j)=av(j)+xx
31 continue
30 continue
c
c.................................... Compute zero-average OOT
s=s/dfloat(n3)
do 6 i=1,n3
v(i)=v(i)-s
6 continue
c
c.................................... Fix transit depth and variance
call robust_aver(njc,y,ave,sig,sigd)
t_ave = dabs(ave)
t_var = sigd*sigd/dfloat(njc)
c
c.................................... Select non-empty bins and calculate
c averages (exclude transit)
nb=0
do 33 j=1,nbq
if(in(j).lt.2) go to 33
if(j.eq.jc) go to 33
nb=nb+1
bv(nb)=av(j)/dfloat(in(j))
33 continue
c
c.................................... Calculate diff_oot
do 34 j=1,nbq
if(in(j).gt.0) av(j)=av(j)/in(j)
34 continue
c
jc1=jc-1
s=0.0d0
k=0
do 35 j=2,jc1
d=av(j)-av(j-1)
s=s+d*d
k=k+1
35 continue
jc2=jc+2
do 36 j=jc2,nbq
d=av(j)-av(j-1)
s=s+d*d
k=k+1
36 continue
oot_diff=s/dfloat(k)
c
c.................................... Calculate variance of the OOT averages
call astat(bv,nb,av_oot,si_oot)
oot_var=si_oot*si_oot
c
c=====================
c Compute DSP
c=====================
c
dsp = t_ave/dsqrt(t_var+oot_var+oot_diff)
if(dsp.gt.99.9d0) dsp=99.9d0
c
c================================
c Compute SNR_ootv, f_ootv
c================================
c
fmin=0.0d0
fmax=30.0d0
nf=3000
df=(fmax-fmin)/dfloat(nf-1)
call fdft(n3,u,v,nf,fmin,df,p,pfr,pam)
call spsnr1(nf,p,fmin,df,frmax1,pmax1,snr1,spd1)
snrt = snr1
fret = frmax1
c
return
end
c
c
subroutine add_flare(n,x)
c---------------------------------------------------------------------
c This routine adds "flares" [positive sequential ouliers] to {x}
c---------------------------------------------------------------------
c PARAMETERS:
c - m :: number of flare events
c - ffac :: amplitude of the flare in the units of xsig
c - idamp :: damping time of the flare in the units of the
c number of data points
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333)
c
m=3
ffac=10.0d0
idamp=10
idamp2=2*idamp
sig_cut=3.0d0
call omit1(sig_cut,x,n,xave,xsig,nclip)
c
k=n/m
k2=k/2
famp=ffac*xsig
fdamp=dfloat(idamp)
do 1 j=1,m
j1=1+(j-1)*k+k2
j2=j1+idamp2
do 2 i=j1,j2
tt=dfloat(i-j1)/fdamp
f=famp*dexp(-tt)
x(i)=x(i)+f
2 continue
1 continue
c
return
end
c
c
subroutine kill_flare(ipr,n,x,nflare)
c---------------------------------------------------------------------
c This routine kills "flares" [positive sequential ouliers] to {x}
c---------------------------------------------------------------------
c NOTES: - Upon return, the number of data points will NOT CHANGE,
c since if a "flare" is found, the corresponding data
c points will be set equal to the average value.
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333)
c
nflare=0
sig_cut=3.0d0
call omit1(sig_cut,x,n,xave,xsig,nclip)
c
rxsig=1.0d0/(sig_cut*xsig)
do 1 i=1,n
d=(x(i)-xave)*rxsig
if(d.lt.1.0d0) go to 1
nflare=nflare+1
x(i)=xave
1 continue
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nflare
if(ipr.eq.1) write(*,11) nflare
if(ipr.eq.2) write(*,12) nflare
10 format('>> nflare [after *PLACE UNKNOWN*] = ',i3)
11 format('>> nflare [after *datain_tfa_four*] = ',i3)
12 format('>> nflare [after *one_rec*] = ',i3)
c
return
end
c
c
subroutine kill_flare_r(ipr,n,x,xf,nflare)
c---------------------------------------------------------------------
c This routine kills "flares" [positive sequential ouliers] in
c *respect to {xf}
c---------------------------------------------------------------------
c NOTES: - Upon return, the number of data points will NOT CHANGE,
c since if a "flare" is found, the corresponding data
c points will be set equal to the corresponding value of
c {xf}.
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),xf(30333),y(30333)
c
nflare=0
sig_cut=3.0d0
do 2 i=1,n
y(i)=x(i)-xf(i)
2 continue
call omit1(sig_cut,y,n,xave,xsig,nclip)
c
rxsig=1.0d0/(sig_cut*xsig)
do 1 i=1,n
d=(y(i)-xave)*rxsig
if(d.lt.1.0d0) go to 1
nflare=nflare+1
x(i)=xf(i)
1 continue
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nflare
if(ipr.eq.1) write(*,11) nflare
if(ipr.eq.2) write(*,12) nflare
10 format('>> nflare [after *PLACE UNKNOWN*] = ',i3)
11 format('>> nflare [after *datain_tfa_four*] = ',i3)
12 format('>> nflare [after *one_rec*] = ',i3)
c
return
end
c
c
subroutine datain_tfa_four(isurvey,nseq,itest,iseed,m_four,
&itemp,fname,sname,ndat,tdat,xdat,nvold,told,xold,
&intemp,avmag,frin0,frin1,pht1,tc0,frootv,amootv,
&phootv,qtr0,qtr1,pht0,igen,std,xsyn,trsyn,
&x_int,camp,we,iflare,ibg,ixy,sig,pfac,ext,
&iwa,fourts,tfats,tfreq0,dep0,tfreq1,dep1,m_opt)
c
c**************************************************************************
c This routine reads in time series and related data and corrects for
c systematics and stellar variability
c**************************************************************************
c
c--------------------------------------------------------------------------
c
c INPUT: - TBD ... or NBD (Never to Be Done) :: At this point (as of
c 2020.09.05 this is just
c too much writing for
c not much advantage ...)
c
c OUTPUT: - {told(i),xold(i); i=1,2,...,nvold} :: ORIGINAL time series
c - {tdat(i); i=1,2,...,ndat}= :: The array {tdat(i)}
c {ttemp(i); i=1,2,...,ntemp} is set equal to the
c TFA template timebase
c - {tdat(i),x_int(i); i=1,2,...,ndat} :: ORIGINAL time series
c interpolated to the
c TFA template timebase
c - {tdat(i),x_tfa(i); i=1,2,...,ndat} :: {x_int}-{fourts}
c
c - {tdat(i),x_four(i); i=1,2,...,ndat} :: {x_int}-{tfats}
c
c - {tdat(i),xdat(i); i=1,2,...,ndat} :: {x_int}-{tfats}-{fourts}
c
c NOTES: -- The Fourier series meant to represent stellar variability,
c by using frequencies 1/T, 2/T, ..., m_four/T [T=total time span]
c NO OTHER frequencies used in this routine.
c -- NO OUTLIER selection in this routine! ALL points are kept!
c -- Time is assumed to be monotonically increasing, i.e.,
c t(i) > t(i-1)
c
c-------------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*3 camp
character*9 sname,idvar
character*25 fname
c
dimension tdat(30333),xdat(30333),t(30333)
dimension told(30333),xold(30333),we(30333)
dimension xsyn(30333),trsyn(30333),x(30333)
dimension fourf(30333),fourts(30333),tfats(30333)
dimension x_int(30333),bgft(30333)
dimension x_p(30333),y_p(30333)
c
common /cc/ m3_tfa,mte0
common /tt/ ttemp(30333),ntemp
common /rc/ m3_rc
c
c************************************************************************
c
idvar=sname
c
c... Set some basic parameters
teps = 1.0d-6
c
c... Read in photometric time series stored in file named *fname*
c
call read_target(itest,iseed,fname,n,t,x,iflare,
&frin0,frin1,tfreq0,dep0,tfreq1,dep1,tc0,
&x_p,y_p,ji,bgft,qtr0,qtr1,pht0,pht1,igen,std,
&xsyn,trsyn,avmag,frootv,amootv,phootv)
c
c... No (X,Y) fit is performed if the input file does not contain
c these time series [safegard against setting ixy=1 in the parameter file]
if(ji.eq.0) ixy=0
c
if(n.eq.0) write(*,22) sname
22 format(' Variable ',a20,' not found !')
if(n.eq.0) return
c
c... Store ORIGINAL INPUT time series
nvold=n
do 15 i=1,nvold
told(i) = t(i)
xold(i) = x(i)
15 continue
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c TFA and TFA+FOUR fits
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
write(*,17)
17 format(/,'Read in TFA templates ... and ...',/,
&'Perform SIMPLE TFA filtering')
c
c Set common timebase, perform classical TFA filtering
c-------------------------------------------------------------
call temp_new(nseq,ji,idvar,itemp,x_p,y_p,bgft,ibg,ixy,
&n,t,x,intemp,x_int)
c
c... Simple TFA filtering
c
call ew_tfa(n,x_int,fourts,tfats,we,sig)
do 16 i=1,n
x(i)=x_int(i)-tfats(i)
16 continue
c
c... Find optimum Fourier order
c
if(m_opt.eq.0) go to 5
write(*,9)
9 format('Search for optimum Fourier order in the RAW data ...')
call four_order_fast(isurvey,ext,n,t,x,m_ord)
m_four=m_ord
write(*,*) 'OPTIMUM Fourier order = ',m_four
5 continue
c
ndat=ntemp
n=ntemp
tmin=ttemp(1)
tmax=ttemp(n)
c
c... Compute sigma-clipped average and standard deviation
c
sig_cut=3.0d0
call omit1(sig_cut,x_int,n,xave,xsig,nclip)
c
tmid=(tmin+tmax)/2.0d0
do 6 i=1,n
t(i) = ttemp(i)-tmid
tdat(i) = t(i)
x_int(i)= x_int(i)-xave
x(i) = x_int(i)
6 continue
c
p0=tmax-tmin
f0=1.0d0/(pfac*p0)
do 8 j=1,m_four
fourf(j)=f0*dfloat(j)
8 continue
c
c............................................... Robust TFA+FOUR filtering
call robust_tfa_four(iwa,n,t,x_int,m_four,fourf,
&fourts,tfats,we,sig)
c
do 1 i=1,n
xdat(i)=x_int(i)-fourts(i)-tfats(i)
1 continue
c
c>>> NOTE: We DO NOT employ single point outlier correction here.
c This is because we may use AR model in the MAIN, that
c could change the outlier status at the edges.
c
call astat(xdat,n,av_x,sig_x)
c
write(*,26) n,sig_x,sig,m3_rc
26 format(/,
&'Number of data points =',i8,/,
&'sig [xdat, no weights, outliers out] =',f8.6,/,
&'sig [xdat, with weights, outliers in] =',f8.6,/,
&'Number of EPD+TFA+FOUR vectors =',i8)
c
c... Reset {tdat}
c
do 59 i=1,n
tdat(i)=tdat(i)+tmid
59 continue
c
if(isurvey.eq.1) go to 10
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c... Save various time series
c
open(54,file='ts_4.dat')
write(54,56) sname,adjustr(camp),avmag,n,m3_tfa,m_four
56 format('#',/,
&'# TARGET: ',a9,/,
&'#',/,'# Campaign: ',a9,/,
&'# : ',f9.5,/,
&'# N_DAT: ',i9,/,
&'# N_TFA: ',i9,/,
&'# N_FOUR: ',i9,/,
&'#',/,
&'#',78('='),/,
&'# TIME',12x,'RAW',8x,'RAW-TFA RAW-FOUR RAW-TFA-FOUR',
&' SYNT',/,
&'#',78('-'))
c
do 58 i=1,n
xx=0.0d0
if(itest.gt.0) xx=xsyn(i)-avmag
write(54,57) tdat(i),x_int(i),x_int(i)-tfats(i),
&x_int(i)-fourts(i),xdat(i),xx
57 format(f14.7,5(2x,f11.7))
58 continue
close(54)
c
10 continue
c
return
end
c
c
subroutine ew_tfa(nt,xa,fourts,tfats,we,sig)
c
c***************************************************************************
c THIS ROUTINE COMPUTES TFA-FILTERED TIME SERIES *
c by using CLASSICAL EQUAL-WEIGHTED Least Squares fit *
c***************************************************************************
c
c Input: nt -- Number of data points of the target time series
c ~~~~~~ ta(i) -- Time values of the target time series
c xa(i) -- Values of the target time series
c
c Output: fourts(i) -- Zero array [outputed only for consistency]
c ~~~~~~~ tfats(i) -- Template function derived from the joint TFA fit
c we(i) -- SUM{we}=nt [just to be consistent with the other
c subroutines]
c sig -- *Unbiased* estimate of the RMS of the residuals
c
c NOTES: -- ZERO POINT value is included in the TFA fit
c -- TFA parameters and arrays are pre-computed, and given
c in the COMMON fields
c
c---------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 idt(995)
c
dimension xa(30333),we(30333)
dimension g(995,996),c(995)
dimension fourts(30333),tfats(30333)
c
common /aa/ xte(995,30333)
common /cc/ m3_tfa,mte0
common /cc1/ idt
common /tt/ ttemp(30333),ntemp
common /rc/ m3_rc
c
c------------------------------------------------------------------------
c
c.................................. Set some parameters, etc.
m3 = m3_tfa
m4 = m3+1
m3_rc= m3
do 2 k=1,nt
fourts(k)=0.0d0
2 continue
c
c.................................. Set weights
rnt=1.0d0
do 1 k=1,nt
we(k)=rnt
1 continue
c
c.................................. Solve equal-weighted LS problem
c
c... Compute TFA normal matrix
c
do 601 i=1,m3
do 602 j=i,m3
s=0.0d0
do 615 k=1,nt
s=s+xte(i,k)*xte(j,k)
615 continue
g(i,j)=s
g(j,i)=s
602 continue
601 continue
c
c... Compute RHS
c
do 303 i=1,m3
s=0.0d0
do 306 k=1,nt
s=s+xa(k)*xte(i,k)
306 continue
g(i,m4)=s
303 continue
c
c... Compute LS solution
c
call gelim1(g,c,m3,IFLAG)
c
c... Compute TFA filter
c
do 26 i=1,nt
s=0.0d0
do 27 j=1,m3
s=s+c(j)*xte(j,i)
27 continue
tfats(i)=s
26 continue
c
c... Compute UNBIASED estimate of the standard deviation
c
s=0.0d0
do 30 k=1,nt
d=xa(k)-tfats(k)
s=s+d*d
30 continue
sig=dsqrt(s/dfloat(nt-m3_tfa))
c
return
end
c
c
subroutine outlier_1g(isurvey,ipr,n,t,x,nout)
c---------------------------------------------------------------------
c This routine selects SINGLE POINT outliers based on 3 neighboring
c points, and then, if an outlier is found, it will be shifted to
c the average of the two neighboring points.
c=====================================================================
c
c NOTES: - Only "middle points" are considered for outlier status
c - It is assumed that the {x} is "nearly" constant with
c zero average.
c - This is a *general* single outlier selection routine.
c It may not work well at sharp variations (e.g., for
c narrow eclipses). If the time series model is known,
c then *outlier_1m* should be used.
c
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),y(30333)
c
sig_cut=3.0d0
c
c... Compute difference time series
c
do 4 i=2,n
i1=i-1
y(i1)=x(i)-x(i1)
4 continue
n1=n-1
c
c... Compute sigma-clipped average and standard deviation
c
call omit1(sig_cut,y,n1,xave,xsig,nclip)
c
if(isurvey.eq.0) open(26,file='outlier.dat')
nout=0
s3=xsig*sig_cut
do 1 i=2,n1
xx=x(i)
i1=i-1
i2=i+1
ave=0.5d0*(x(i1)+x(i2))
d=dabs(xx-ave)
if(d.lt.s3) go to 2
if(dabs(x(i2)-x(i1)).gt.d) go to 2
nout=nout+1
x(i)=ave
if(isurvey.eq.0) write(26,3) t(i),xx,x(i)
3 format(f12.5,2(2x,f10.7))
2 continue
1 continue
if(isurvey.eq.0) close(26)
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nout
if(ipr.eq.1) write(*,11) nout
if(ipr.eq.2) write(*,12) nout
10 format('>> nout [after *PLACE UNKNOWN*] = ',i3)
11 format('>> nout [after *datain_tfa_four*] = ',i3)
12 format('>> nout [after *one_rec*] = ',i3)
c
return
end
c
c
subroutine outlier_1m(isurvey,ipr,n,t,x,fx,nout)
c---------------------------------------------------------------------
c This routine selects SINGLE POINT outliers based on 3 neighboring
c points. The outlier status is controlled by the local difference
c between the data {x} and the model {fx}
c=====================================================================
c
c NOTES: - Only "middle points" are considered for outlier status
c
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),fx(30333),y(30333)
c
sig_cut=3.0d0
cut=3.0d0
n1=n-1
c
c... Compute difference time series
c
s1=0.0d0
do 4 i=1,n
y(i)=x(i)-fx(i)
s1=s1+y(i)
4 continue
s1=s1/dfloat(n)
c
c... Compute sigma-clipped average and standard deviation of {y}
c
call omit1(sig_cut,y,n1,xave,xsig,nclip)
c
s=0.0d0
do 14 i=1,n
y(i)=y(i)-s1
s=s+y(i)
14 continue
c
if(isurvey.eq.0) open(26,file='outlier.dat')
nout=0
s3=xsig*cut
do 1 i=2,n1
xx=x(i)
i1=i-1
i2=i+1
c
c........................................ Neighbors CANNOT be outliers!
if(dabs(y(i)).lt.s3) go to 2
if(dabs(y(i1)).ge.s3) go to 2
if(dabs(y(i2)).ge.s3) go to 2
c........................................
nout=nout+1
if(isurvey.eq.0) write(26,3) t(i),x(i),fx(i)
x(i)=fx(i)
3 format(f12.5,2(2x,f10.7))
2 continue
1 continue
if(isurvey.eq.0) close(26)
c
if(ipr.lt.0) return
if(ipr.eq.0) write(*,10) nout
if(ipr.eq.1) write(*,11) nout
if(ipr.eq.2) write(*,12) nout
10 format('>> nout [after *PLACE UNKNOWN*] = ',i3)
11 format('>> nout [after *datain_tfa_four*] = ',i3)
12 format('>> nout [after *one_rec*] = ',i3)
c
c
return
end
c
c
subroutine bls_white(n,t,x,nf1,fmin1,df1,nb,qmi,qma,
&nsing,p0_sing,sname,nb_sp,ndsp,nbb,ntw,fw,dw,snrw)
c---------------------------------------------------------------------
c This routine performs successive BLS prewhitening
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*8 spn,lcn,pan
character*9 sname
c
dimension t(30333),x(30333),y(30333)
dimension p(155555),xbt(30333)
dimension u(30333),v(30333)
dimension fw(100),dw(100),snrw(100)
dimension intr(30333)
c
c=================================================================================
c
write(*,*) '>>>>> SUCCESSIVE BLS signal search in *bls_white*'
c
open(25,file='multi_bls.dat')
write(25,3)
write(*,3)
3 format('#',67('='),/,
&'# iw f[1/d] depth qtran Tc N_intr ',
&'SNR N_event',/,
&'#',67('-'))
c
call astat(x,n,aver,disper)
do 6 i=1,n
y(i)=x(i)-aver
6 continue
c
do 1 iw=1,ntw
c
iw1=iw-1
c
c................................................................. SP, LC, TR parameter file names
call sp_lc_name(iw1,spn,lcn,pan)
c
c................................................................. BLS
call eebls(n,t,y,u,v,nf1,fmin1,df1,nb,qmi,qma,
&p,bper,bpow,depth,qtran,in1,in2)
c
c................................................................. SP polfit
call sppolfit_w(ndsp,nf1,p)
c
c................................................................. Saving SP
call save_sp_iw(nb_sp,nf1,fmin1,df1,p,spn)
c
c................................................................. SNR, frmax0
call spsnr1(nf1,p,fmin1,df1,frmax0,pmax0,sde0,spd0)
c
c................................................................. SYNTHETIC signal
c
c....................... Accurate period
p0=1.0d0/frmax0
call ac_per(p0,df1,n,t,y,nbb,qmi,qma,pac,sig)
p0=pac
c
c....................... Use SINGLE transit period if warranted
if(nsing.eq.1.and.iw.eq.1) p0=p0_sing
c
c....................... Transit parameters
call foldts4_new(p0,qmi,qma,nbb,n,t,y,xbt,
&intr,tdepth,ring,qtran,epc,nintr,sig4)
c
c................................................................. Saving LC
call save_lc_iw(n,t,y,xbt,lcn)
c
c................................................................. Saving folding parameters
call save_pa_iw(sname,p0,epc,pan)
c
c................................................................. WHITENING
s=0.0d0
do 4 i=1,n
y(i)=y(i)-xbt(i)
s=s+y(i)
4 continue
s=s/dfloat(n)
do 5 i=1,n
y(i)=y(i)-s
5 continue
c
fw(iw) = 1.0d0/p0
dw(iw) = tdepth
snrw(iw) = sde0
c
f0=fw(iw)
tcen=epc
call events(n,t,f0,tcen,qtran,nev,ndit,itoc)
c
write(*,2) iw,fw(iw),tdepth,qtran,epc,nintr,sde0,nev
write(25,2) iw,fw(iw),tdepth,qtran,epc,nintr,sde0,nev
2 format(2x,i2,1x,f8.5,2x,f9.6,2x,f7.5,2x,f13.7,2x,i4,
&2x,f5.1,1x,i3)
c
1 continue
write(*,7)
write(25,7)
7 format('#',67('-'),/,'#')
close(25)
c
return
end
c
c
subroutine sp_lc_name(iw,spn,lcn,pan)
c---------------------------------------------------------------------
c Create sequential file names as given by *iw* :
c
c spn = Spectrum from the iw-th prewhitening step
c lcn = Light curve from the iw-th prewhitening step
c pan = Folding parameter file appropriate for the above datasets
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*2 ctemp
character*8 spn,lcn,pan
c
iw0=iw
if(iw.gt.99) write(*,1)
1 format('File number cannot be greater than 99 !')
if(iw.gt.99) iw=99
c
write(ctemp,'(i2)') iw
if(iw.lt.10) spn='sp0'//trim(adjustl(ctemp))//'.dat'
if(iw.ge.10) spn='sp'//trim(adjustl(ctemp))//'.dat'
if(iw.lt.10) lcn='lc0'//trim(adjustl(ctemp))//'.dat'
if(iw.ge.10) lcn='lc'//trim(adjustl(ctemp))//'.dat'
if(iw.lt.10) pan='pa0'//trim(adjustl(ctemp))//'.dat'
if(iw.ge.10) pan='pa'//trim(adjustl(ctemp))//'.dat'
iw=iw0
c
return
end
c
c
subroutine one_rec(isurvey,iew,jsm,p0,nt,ta,xi,xa,nbb,
&qmi1,qma1,m_four,iso,intr,nintr,tdepth,ring,rtran,tcen,
&xft,xtr,we,sig,nout,tfa_r,four_r)
c
c***************************************************************************
c SINGLE PERIODIC TRANSIT SIGNAL RECONSTRUCTION *
c***************************************************************************
c
c Input: iew -- 1/0/2 use/do not use equal weights and exact zero
c ~~~~~~ weights (2)
c jsm -- 0 = no reconstruction; 1 = reconstruction
c p0 -- Period of the transit signal
c nt -- Number of data points
c {ta} -- Timebase ( {ta}={ttemp} )
c {xi} -- Original time series, interpolated to the TFA
c template timebase ( {xi}={x_int} )
c {xa} -- Some systematics-filtered version of the target
c time series [to compute initial estimate on the
c transit signal {xft}; {xa}={x}]
c {we} -- Initial weights
c nbb -- Number of bins for "eebls"
c qmi1 -- Minimum q_tran for "eebls"
c qma1 -- Maximum q_tran for "eebls"
c m_four-- Fourier order
c iso -- 0/1 - no/yes single outlier correction
c
c Output: {intr} -- Pointer array: 1 if 'intransit', 0 if 'oot'
c ~~~~~~~ nintr -- INTRANSIT (between t1 and t4) data point number
c tdepth -- Transit depth
c ring -- t12/t14
c rtran -- t14/P
c tcen -- T_c [time at the center of the transit]
c {xft} -- Best-fitting transit signal
c {xtr} -- {xi}-{tfats}-{fourts}
c [{xft} is the best fit transit signal to {xtr}]
c {we} -- Weights from the best fitting full model
c sig -- **BIASED** estimate of the standard deviation
c of the residuals of the time series with the
c TFA+FOUR+TRANSIT components subtracted from the
c input time series.
c nout -- Number of single outliers
c {tfa_r} -- TFA vector [including also the EPD parameters]
c resulting from the FULL model fit
c {four_r} -- FOUR vector resulting from the FULL model fit
c
c
c NOTES: -- {xi} = {x_int} in MAIN
c -- {xa} = {x} in MAIN ({xdat} in "datain_tfa_four")
c [i.e., {xa} is the "non-reconstructed" time
c series]
c -- TFA parameters and arrays are pre-computed, and given
c in the COMMON fields
c -- The full reconstruction includes:
c transit + systematics + stellar variability;
c {xte} includes systematics + stellar variability
c -- There are altogether "m3_rc" correcting vectors.
c We add the updated approximation of the transit signal
c {xft} at each iteration to the above set of correcting
c vectors. Thereby we enable the full reconstruction.
c -- Since during each iteration step the *original* data
c are fitted, the single outliers reported are NOT
c CUMMULATIVE.
c -- Iteration on the complete solution starts with the
c weights {we} used in the pre-BLS phase.
c -- A transit solution is searched for, even if the best
c solution results in a positive dip.
c
c------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension xi(30333),ta(30333),xa(30333),we(30333)
dimension xft(30333),xtr(30333),xw(30333),si(30333)
dimension wa(30333),tfa_r(30333),four_r(30333)
dimension g(995,996),c(995)
dimension intr(30333)
c
common /aa/ xte(995,30333)
common /cc/ m3_tfa,mte0
common /rc/ m3_rc
c
c------------------------------------------------------------------------
c
c.................................. Set parameters for the densely-sampled fit
n_ph = 200
r_q = 1.5d0
c
c.................................. Set iteration parameters, etc.
s3 = 1.0d0
nit = 15
ita = 0
nit5 = 5
if(jsm.eq.1.or.iso.eq.0) nit5 = 1
eps = 0.001d0
sig0 = 1.0d33
m3 = m3_rc+1
m4 = m3+1
m5 = m3_rc
c--------------------------
c
c... First estimate of the transit signal from {xa}
c
do 102 i=1,nt
xtr(i)=xa(i)
wa(i)=1.0d0
102 continue
c
c... Iterate on the NON-RECONSTRUCTED time series for single outliers
c
do 200 it=1,nit5
call foldts4_new(p0,qmi1,qma1,nbb,nt,ta,xtr,xft,
&intr,tdepth,ring,rtran,epc,nintr,sig4)
if(iso.eq.1) call outlier_1m(isurvey,-1,nt,ta,xtr,xft,nout)
200 continue
tcen=epc
sig=sig4
c
if(jsm.eq.0) return
if(iew.eq.2) call slide_rms(m_four,nt,ta,xtr,si,wa)
c
c.................................. ITERATION on {xft}
c
write(*,101)
101 format('>>>>> ITERATION in *one_rec*')
c
do 100 it=1,nit
c
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 100
ita=ita+1
sig0=sig
c
c... Add estimated transit signal to the TFA template set
c
do 32 k=1,nt
xte(m3,k)=xft(k)
32 continue
c
c... Compute TFA normal matrix
c
do 601 i=1,m3
do 602 j=i,m3
s=0.0d0
do 615 k=1,nt
s=s+we(k)*xte(i,k)*xte(j,k)
615 continue
g(i,j)=s
g(j,i)=s
602 continue
601 continue
c
c... Compute RHS
c
do 303 i=1,m3
s=0.0d0
do 306 k=1,nt
s=s+we(k)*xi(k)*xte(i,k)
306 continue
g(i,m4)=s
303 continue
c
c... Compute LS solution
c
call gelim1(g,c,m3,IFLAG)
c
c... Compute transit signal
c {xtr} = {xi}-{TFA}-{FOUR}
c and standard deviation of the residual
c {res} = {xi}-{TFA}-{FOUR}-{TRANSIT}
c
d=0.0d0
do 26 k=1,nt
s=0.0d0
do 27 j=1,m5
s=s+c(j)*xte(j,k)
27 continue
xx=xi(k)-s
yy=xx-c(m3)*xte(m3,k)
d=d+yy*yy
xtr(k)=xx
xw(k)=yy
26 continue
c
c... Compute BIASED estimate of the equal-weighted standard deviation
c
sig=dsqrt(d/dfloat(nt))
if(iew.eq.1) go to 99
c
c.................................. UPDATE {we}
cc=s3*sig
c2=cc*cc
s=0.0d0
do 5 k=1,nt
dd=xw(k)
d2=dd*dd
d2=1.0d0/(c2+d2)
d2=wa(k)*d2
we(k)=d2
s=s+d2
5 continue
s=1.0d0/s
do 6 k=1,nt
we(k)=s*we(k)
6 continue
c
99 continue
c
c.................................. UPDATE {xft}
c
call foldts4_w(p0,qmi1,qma1,nbb,nt,ta,xtr,xft,
&we,intr,tdepth,ring,rtran,epc,nintr,sig4)
tcen=epc
sig=sig4
if(iew.eq.2) call slide_rms(m_four,nt,ta,xtr,si,wa)
c
100 continue
c
write(*,1) ita
1 format('Actual iteration number in *one_rec* :',i3)
c
c.................................. Compute FOUR and TFA components
c.................................. 1 - n0 : TFA templates
c n2 - n3 : FOURIER
m2 = 2*m_four
n0 = m3_tfa
n2 = n0+1
n3 = n0+m2
c
c... Compute TFA filter
c
do 10 i=1,nt
s=0.0d0
do 11 j=1,n0
s=s+c(j)*xte(j,i)
11 continue
tfa_r(i)=s
10 continue
c
c... Compute best fitting FOUR time series
c
do 12 i=1,nt
s=0.0d0
do 13 j=n2,n3
s=s+c(j)*xte(j,i)
13 continue
four_r(i)=s
12 continue
c
if(iso.eq.0) go to 40
call outlier_1m(isurvey,2,nt,ta,xtr,xft,nout)
call foldts4_w(p0,qmi1,qma1,nbb,nt,ta,xtr,xft,
&we,intr,tdepth,ring,rtran,epc,nintr,sig4)
tcen=epc
sig=sig4
c
40 continue
c
c... Generate and save densely sampled transit fit
c
ktr=1
d00=tdepth
q00=rtran
r00=ring
call save_trc(ktr,d00,q00,r00)
c
return
end
c
c
subroutine phas05(time,e0,f0,ph)
c-------------------------------------------------------------------
c This routine computes phase in [-0.5,0.5]
c-------------------------------------------------------------------
implicit real*8 (a-h, o-z)
implicit integer*4 (i-n)
c
ph = f0*(time-e0)
ph = ph - idint(ph)
if(ph.gt. 0.5d0) ph=ph-1.0d0
if(ph.lt.-0.5d0) ph=ph+1.0d0
c
return
end
c
c
subroutine trapu(n,t,p0,e0,q,r,tru)
c===================================================================
c This routine generates periodic transit time series on {t}
c-------------------------------------------------------------------
c
c DATE: 2019.12.23
c
c INPUT: - n :: Number of data points
c - {t} :: Time values ["n" time values altogether]
c - p0 :: Period
c - e0 :: Epoch of the center of the transit
c - q :: t14/p0
c - r :: t12/t14
c
c OUTPUT: - {tru} :: "n" transit function values generated on the
c timebase as given by {t}
c
c NOTES: - Flat OOT and U-shaped bottom and linear ingress/egress
c - OOT = 0.0
c - Center of the transit has a depth of: -1.0
c - Parameter "r" is defined relative to the total transit
c duration and NOT to the period!
c
c-------------------------------------------------------------------
c
implicit real*8 (a-h, o-z)
implicit integer*4 (i-n)
c
dimension t(30333),tru(30333)
c
c... Set basic transit parameters
c
cu = 0.10d0
eps = 1.0d-10
c
f0=1.0d0/p0
cu1= cu-1.0d0
rq = r*q
ph1=-q/2.0d0
ph4=ph1+q
ph2=ph1+rq
ph3=ph4-rq
xfac=2.0d0/(ph3-ph2)
c
if(rq.gt.eps) go to 1
c
c--------------------------------
c Infinite ingress slope
c--------------------------------
do 2 k=1,n
c
time=t(k)
call phas05(time,e0,f0,ph)
tr = 0.0d0
if(ph.lt.ph1) go to 9
if(ph.gt.ph4) go to 9
xx=ph*xfac
x2=xx*xx
tr = -1.0d0+cu*(1.0d0-dsqrt(1.0d0-x2))
9 continue
tru(k)=tr
2 continue
return
c
1 continue
c--------------------------------
c Finite ingress slope
c--------------------------------
r12=cu1/rq
c
do 11 k=1,n
c
time=t(k)
call phas05(time,e0,f0,ph)
tr = 0.0d0
if(ph.lt.ph1) go to 12
if(ph.gt.ph4) go to 12
if(ph.lt.ph2) go to 13
if(ph.gt.ph3) go to 14
c
c..................................... Bottom of the transit ph2 < ph < ph3
xx=ph*xfac
x2=xx*xx
tr = -1.0d0+cu*(1.0d0-dsqrt(1.0d0-x2))
go to 12
c
c..................................... Ingress/Egress parts
13 continue
tr = r12*(ph-ph1)
go to 12
14 continue
tr = cu1-r12*(ph-ph3)
c
12 continue
tru(k)=tr
11 continue
c
return
end
c
c
subroutine foldts4_new(p0,qmi1,qma1,nb,n,t,x,tru,
&intr,tdepth,ring,rtran,epc,nintr,sig4)
c
c**********************************************************************
c THIS ROUTINE FITS an U-shaped trapezoidal signal with FLAT OOTV *
c to a time series *
c**********************************************************************
c
c INPUT: p0 -- folding period
c qmi1 -- Minimum fractional transit length used in the
c trapezoidal fit
c qma1 -- Maximum fractional transit length used in the
c trapezoidal fit
c nb -- Number of bins in the BLS approximation
c n -- Number of data points of the input time series
c {t} -- Time values of the time series
c {x} -- Data values of the time series
c
c OUTPUT: {tru} -- Best transit fit
c {intr} -- Pointer array: 1 if 'intransit', 0 if
c 'out-of-transit'
c tdepth -- Depth of the transit (deepest value of the
c transit}
c ring -- t12/P (=t34/P)
c rtran -- t14/P
c epc -- Epoch of the center of the transit
c nintr -- INTRANSIT (between t1 and t4) data point number
c sig4 -- **BIASED** estimate of the RMS of the residuals
c
c NOTES: - We assume that the moments of time are monotonic, i.e.,
c t(i) > t(i-1)
c
c------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),tru(30333),tr0(30333)
dimension u(30333),v(30333),p(155555)
dimension intr(30333)
c
c------------------------------------------------------------------------
c
c... Wired in parameters [ONLY FOR THIS ROUTINE]
c
exq = 0.3d0
nc = 7
ntrl = 7
ning = 7
eps = 1.0d-5
eps0 = -1.0d-20
itm = 50
qfac = 0.321d0
c
c exq = Total transit duration (t4-t1) is scanned in the range
c of [qtran*(1-exq),qtran*(1+exq)], where *qtran* is
c the estimate given by the routine 'eebls'
c .......... NOTE: we use the SAME factor for the EPOCH
c nc = Center of transit scan number
c ntrl = Transit length scan number
c ning = Ingress length scan number
c eps = When the step size of the central transit parameter
c becomes smaller than *eps*, iteration stops. "eps" is
c measured in [d]
c eps0 = OOT level in making distinction between the transit and
c OOT of the synthetic transit signal as generated by "trapu"
c [the output signal of that routine has a zero OOT level]
c itm = Maximum number of iterations
c qfac = Scanning regime limiting parameter. It should be in (0,1).
c Smaller values of it lead to the narrower parameter regimes
c to be refined in the next iteration.
c
c------------------------------------------------------------------------
c
c-----------------
c... Set some parameters just to avoid
c 'might be used uninitialized' messages during compilation
c
c0b=0.0d0
c1b=1.0d0
rb=0.0d0
qb=0.03d0
eb=0.0d0
c-----------------
c
c... Find transit parameters in the BLS approximation
c
nf1 = 1
fmin1 = 1.0/p0
df1 = 0.0d0
c
call eebls(n,t,x,u,v,nf1,fmin1,df1,nb,qmi1,qma1,
&p,bper,bpow,depth,qtran,in1,in2)
c
c... Set EPOCH [ here we assume that "eebls" uses {t(i)-tmin} ]
c
tmin=t(1)
tmax=t(n)
c
if(in1.gt.in2) in2=in2+nb
phmid=dfloat(in1+in2)/dfloat(2*nb)
tmid=p0*phmid
e=tmin+tmid
c
c****************************
c Start scanning
c****************************
c
er = qtran*p0*exq
e1 = e - er
e2 = e + er
de = (e2-e1)/dfloat(nc-1)
qr = qtran*exq
q1 = qtran - qr
q2 = qtran + qr
dq = (q2-q1)/dfloat(ntrl-1)
rr = 0.25d0
r1 = 0.25d0-rr
r2 = 0.25d0+rr
dr = (r2-r1)/dfloat(ning-1)
sdm = 1.0d30
it = 0
c
30 continue
it = it +1
c
c... Center of transit - e0
c
do 20 ic=1,nc
c
e = e1 + de*dfloat(ic-1)
c
c... Total transit length - q
c
do 21 iq=1,ntrl
c
q = q1 + dq*dfloat(iq-1)
c
c... Ingress length - r
c
do 22 iing=1,ning
c
r = r1 + dr*dfloat(iing-1)
c
call trapu(n,t,p0,e,q,r,tru)
call tru_fit(n,tru,x,c0,c1,sdev)
c
if(sdev.gt.sdm) go to 22
sdm = sdev
eb = e
qb = q
rb = r
c0b = c0
c1b = c1
c
22 continue
21 continue
20 continue
c
c... Check convergence criterion
c
if(dabs(eb-e).lt.eps) go to 31
e=eb
c
c... Iteration has NOT converged: update search regime
c
if(it.gt.itm) write(*,23) sdm
23 format('Iteration in *foldts4_new* has not converged!',/,
&'Sigma=',f11.5,'Computation is continued anyway ...')
if(it.gt.itm) go to 31
c
er = er*qfac
e1 = eb - er
e2 = eb + er
de = (e2-e1)/dfloat(nc-1)
qr = qr*qfac
q1 = qb - qr
if(q1.lt.0.0d0) q1=0.0d0
q2 = qb + qr
dq = (q2-q1)/dfloat(ntrl-1)
rr = rr*qfac
r1 = rb - rr
if(r1.lt.0.0d0) r1=0.0d0
r2 = rb +rr
if(r2.gt.0.5d0) r2=0.5d0
dr = (r2-r1)/dfloat(ning-1)
c
go to 30
c
31 continue
c
c... Iteration has converged:
c - assign output parameter names
c - compute best transit function
c
call trapu(n,t,p0,eb,qb,rb,tru)
call tru_fit(n,tru,x,c0b,c1b,sdev)
sig4=sdev
d1=1.0d30
d2=-d1
s=0.0d0
do 36 i=1,n
tr_u=tru(i)
tr0(i)=tr_u
tr_u=c0b+c1b*tr_u
tru(i)=tr_u
d1=dmin1(d1,tr_u)
d2=dmax1(d2,tr_u)
36 continue
tdepth = d2-d1
if(c1b.lt.0.0d0) tdepth=-tdepth
ring = rb
rtran = qb
epc = eb
c
c... Compute intransit data point number and intransit pointer array {intr}
c
nintr = 0
do 35 i=1,n
intr(i)=0
if(tr0(i).gt.eps0) go to 35
nintr = nintr + 1
intr(i)=1
35 continue
c
return
end
c
c
subroutine foldts4_w(p0,qmi1,qma1,nb,n,t,x,tru,
&we,intr,tdepth,ring,rtran,epc,nintr,sig4)
c
c**********************************************************************
c THIS ROUTINE FITS an U-shaped trapezoidal signal with FLAT OOTV *
c to a time series. WEIGHTED Least Squares is used *
c**********************************************************************
c
c INPUT: p0 -- folding period
c qmi1 -- Minimum fractional transit length used in the
c trapezoidal fit
c qma1 -- Maximum fractional transit length used in the
c trapezoidal fit
c nb -- Number of bins in the BLS approximation
c n -- Number of data points of the input time series
c {t} -- Time values of the time series
c {x} -- Data values of the time series
c {we} -- Point-by-point weight
c
c OUTPUT: {tru} -- Best transit fit
c {intr} -- Pointer array: 1 if 'intransit', 0 if
c 'out-of-transit'
c tdepth -- Depth of the transit (deepest value of the
c transit}
c ring -- t12/P (=t34/P)
c rtran -- t14/P
c epc -- Epoch of the center of the transit
c nintr -- INTRANSIT (between t1 and t4) data point number
c sig4 -- **BIASED** estimate of the RMS of the residuals
c
c NOTES: - We assume that the moments of time are monotonic, i.e.,
c t(i) > t(i-1)
c
c------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),tru(30333),tr0(30333)
dimension u(30333),v(30333),we(30333),p(155555)
dimension intr(30333)
c
c------------------------------------------------------------------------
c
c... Wired in parameters [ONLY FOR THIS ROUTINE]
c
exq = 0.3d0
nc = 7
ntrl = 7
ning = 7
eps = 1.0d-5
eps0 = -1.0d-20
itm = 50
qfac = 0.321d0
c
c exq = Total transit duration (t4-t1) is scanned in the range
c of [qtran*(1-exq),qtran*(1+exq)], where *qtran* is
c the estimate given by the routine 'eebls'
c .......... NOTE: we use the SAME factor for the EPOCH
c nc = Center of transit scan number
c ntrl = Transit length scan number
c ning = Ingress length scan number
c eps = When the step size of the central transit parameter
c becomes smaller than *eps*, iteration stops. "eps" is
c measured in [d]
c eps0 = OOT level in making distinction between the transit and
c OOT of the synthetic transit signal as generated by "trapu"
c [the output signal of that routine has a zero OOT level]
c itm = Maximum number of iterations
c qfac = Scanning regime limiting parameter. It should be in (0,1).
c Smaller values of it lead to the narrower parameter regimes
c to be refined in the next iteration.
c
c------------------------------------------------------------------------
c
c-----------------
c... Set some parameters just to avoid
c 'might be used uninitialized' messages during compilation
c
c0b=0.0d0
c1b=1.0d0
c-----------------
c
c... Find transit parameters in the BLS approximation
c
nf1 = 1
fmin1 = 1.0/p0
df1 = 0.0d0
c
call weebls_new(n,t,x,we,u,v,nf1,fmin1,df1,nb,qmi1,qma1,
&p,bper,bpow,depth,qtran,in1,in2)
c
c... Set EPOCH [ here we assume that "eebls" uses {t(i)-tmin} ]
c
tmin=t(1)
tmax=t(n)
c
if(in1.gt.in2) in2=in2+nb
phmid=dfloat(in1+in2)/dfloat(2*nb)
tmid=p0*phmid
e=tmin+tmid
c
c****************************
c Start scanning
c****************************
c
er = qtran*p0*exq
e1 = e - er
e2 = e + er
de = (e2-e1)/dfloat(nc-1)
qr = qtran*exq
q1 = qtran - qr
q2 = qtran + qr
dq = (q2-q1)/dfloat(ntrl-1)
rr = 0.25d0
r1 = 0.25d0-rr
r2 = 0.25d0+rr
dr = (r2-r1)/dfloat(ning-1)
sdm = 1.0d30
it = 0
c
30 continue
it = it +1
c
c... Center of transit - e0
c
do 20 ic=1,nc
c
e = e1 + de*dfloat(ic-1)
c
c... Total transit length - q
c
do 21 iq=1,ntrl
c
q = q1 + dq*dfloat(iq-1)
c
c... Ingress length - r
c
do 22 iing=1,ning
c
r = r1 + dr*dfloat(iing-1)
c
call trapu(n,t,p0,e,q,r,tru)
call tru_fit_w(n,we,tru,x,c0,c1,sdev)
c
if(sdev.gt.sdm) go to 22
sdm = sdev
eb = e
qb = q
rb = r
c0b = c0
c1b = c1
c
22 continue
21 continue
20 continue
c
c... Check convergence criterion
c
if(dabs(eb-e).lt.eps) go to 31
e=eb
c
c... Iteration has NOT converged: update search regime
c
if(it.gt.itm) write(*,23) sdm
23 format('Iteration in *foldts4_new* has not converged!',/,
&'Sigma=',f11.5,'Computation is continued anyway ...')
if(it.gt.itm) go to 31
c
er = er*qfac
e1 = eb - er
e2 = eb + er
de = (e2-e1)/dfloat(nc-1)
qr = qr*qfac
q1 = qb - qr
if(q1.lt.0.0d0) q1=0.0d0
q2 = qb + qr
dq = (q2-q1)/dfloat(ntrl-1)
rr = rr*qfac
r1 = rb - rr
if(r1.lt.0.0d0) r1=0.0d0
r2 = rb +rr
if(r2.gt.0.5d0) r2=0.5d0
dr = (r2-r1)/dfloat(ning-1)
c
go to 30
c
31 continue
c
c... Iteration has converged:
c - assign output parameter names
c - compute best transit function
c
call trapu(n,t,p0,eb,qb,rb,tru)
call tru_fit_w(n,we,tru,x,c0b,c1b,sdev)
sig4=sdev
d1=1.0d30
d2=-d1
s=0.0d0
do 36 i=1,n
tr_u=tru(i)
tr0(i)=tr_u
tr_u=c0b+c1b*tr_u
tru(i)=tr_u
d1=dmin1(d1,tr_u)
d2=dmax1(d2,tr_u)
36 continue
tdepth = d2-d1
if(c1b.lt.0.0d0) tdepth=-tdepth
ring = rb
rtran = qb
epc = eb
c
c... Compute intransit data point number and intransit pointer array {intr}
c
nintr = 0
do 35 i=1,n
intr(i)=0
if(tr0(i).gt.eps0) go to 35
nintr = nintr + 1
intr(i)=1
35 continue
c
return
end
c
c
subroutine ac_per(p0,df1,n,t,x,nb,qmi1,qma1,pac,sig)
c-------------------------------------------------------------
c This routine computes ACCURATE TRANSIT PERIOD
c-------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),tru(30333)
dimension intr(30333)
c
nit = 50
eps = 1.0d-4
nf = 10
tot = t(n)-t(1)
c
fmin_min = 1.001d0/tot
edf = eps/tot
f0 = 1.0d0/p0
facf = 1.54321d0
fmin = f0-facf*df1
if(fmin.lt.fmin_min) fmin=fmin_min
fmax = fmin+2.0d0*facf*df1
it = 0
f1 = f0
f2 = f1
sigm = 1.0d30
c
1 continue
it=it+1
if(it.gt.nit) go to 2
df=(fmax-fmin)/dfloat(nf-1)
c
do 4 j=1,nf
ff=fmin+df*dfloat(j-1)
pp=1.0d0/ff
call foldts4_new(pp,qmi1,qma1,nb,n,t,x,tru,
&intr,tdepth,ring,rtran,epc,nintr,sig4)
if(sig4.gt.sigm) go to 4
sigm=sig4
f2=ff
4 continue
if(dabs(f2-f1).lt.edf) go to 2
fmin=f2-facf*df
if(fmin.lt.fmin_min) fmin=fmin_min
fmax=fmin+2.0d0*facf*df
f1=f2
go to 1
2 continue
pac=1.0d0/f2
sig=sigm
c
return
end
c
c
subroutine slide_rms(m,n,t,x,si,wa)
c
c*************************************************
c This one is cumbersome, because of the continuous ZERO weight
c region search. The aim is to limit the size of the ZERO weight
c intervals. The effect of this constraint on the final result
c is unclear. THEREFORE, we use the simpler ZERO WEIGHT routine
c (without the "old" extension) for the time being.
c
c DATE: 2020.02.26
c
c*************************************************
c---------------------------------------------------------------------
c This is to compute binary (0 or 1) weight function {wa}
c---------------------------------------------------------------------
c
c INPUT: m - Order of the Fourier series
c n - Number of data points
c {t} - Moments of time (n values)
c {x} - Time series values
c
c OUTPUT: {si} - Normalized standard deviations within the
c sliding boxes (n values)
c {wa} - Binary weights (0 or 1 for large/small standard
c deviations {si}; n values)
c
c NOTES: - This routine computes sliding RMS on a time series
c - Assign binary weights to the high RMS intervals
c - We assume nearly equidistant distribution in time
c fraction of the total number of data points.
c - We assume monotonically increasing moments of time
c - Endpoint weights (wa(1) and wa(n)) are set equal to 1
c - We assume that the moments of time are monotonic.
c
c=====================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),wa(30333)
dimension v(30333),si(30333)
dimension is(30333),ie(30333)
c
c............................... Set some basic parameters
asig = 2.0d0
wcut = 2.0d0
dt_max = dmin1(1.50d0,0.75d0*(t(n)-t(1))/dfloat(m))
dt = 0.50d0
eps = 1.0d-6
c
c asig - sigma-clipping parameter for "omit1" to compute
c outlier-free value of the average standard deviation
c of the standard deviations within the sliding box
c
c wcut - If the given item of data has an associated sliding box
c standard deviation greater than wcut*xave0, then the
c associated weight array element (wa(k)) is set to ZERO.
c
c dt_max - Maximum size (in [days]) of a single time interval
c assigned with ZERO weight. We take 3/4 of the shortest
c period fitted, to avoid overshooting in the {wa} windowed
c intervals.
c
c dt - Sliding box size [days]
c
c eps - Values in {wa} are considered jumping if
c |wa(i)-wa(i-1)| > eps
c
c............................... Set index limits for end points
c [because of the finite width of the
c sliding box, the first RMS can be
c computed only in the middle of the box]
dt2=dt/2.0d0
n1=0
tmin=t(1)
tmax=t(n)
5 continue
n1=n1+1
if(t(n1)-tmin.lt.dt2) go to 5
n2=n
6 continue
n2=n2-1
if(tmax-t(n2).lt.dt2) go to 6
nw=(n1+n-n2)/2
n3=2*nw+1
c
c................................... Scan {t,x} with a sliding box of a
c length {2*nw}
do 1 i=n1,n2
i1=i-nw
i2=i+nw
do 2 j=i1,i2
k=j-i1+1
v(k)=x(j)
2 continue
call astat(v,n3,xave,xsig)
si(i)=xsig
1 continue
c................................... End sections are set constant
c (corresponding to the first/last
c computed values of {si})
si1=si(n1)
do 50 i=1,n1
si(i)=si1
50 continue
n4=n2+1
si2=si(n2)
do 51 i=n4,n
si(i)=si2
51 continue
c................................... Compute robust average standard deviations
c derived from the sliding box analysis
call omit1(asig,si,n,xave0,xsig0,nclip0)
fac=1.0d0/xave0
c................................... Normalize {si} with the above average
c standard deviation
do 52 i=1,n
si(i)=fac*si(i)
52 continue
c
c................................... Compute bimodal weights {wa}
n4=n-1
wa(1)=1.0d0
wa(n)=1.0d0
do 10 i=2,n4
wa(i)=1.0d0
if(si(i).gt.wcut) wa(i)=0.0d0
10 continue
c
c*******************************************************
c Find array indices of ZERO {wa} intervals
c - {is} : starting indices of the zero intervals
c - {ie} : ending indices of the zero intervals
c*******************************************************
c
w1=wa(1)
ns=0
ne=0
do 11 i=2,n
w2=wa(i)
if(dabs(w2-w1).lt.eps) go to 12
if(w2-w1.lt.-0.5d0) go to 13
if(w2-w1.gt.0.5d0) go to 14
13 continue
ns=ns+1
is(ns)=i
go to 111
14 continue
ne=ne+1
ie(ne)=i-1
111 continue
12 continue
w1=w2
11 continue
if(ns.ne.ne) write(*,*) 'ns MUST be equal to ne !!'
if(ns.ne.ne) ns=ne
c
c........................................ Shrink high-sigma region to dt_max
do 15 j=1,ns
tt=t(ie(j))-t(is(j))
if(tt.lt.dt_max) go to 15
rat=dt_max/tt
n_mid=1+idint(0.5d0*(dfloat(is(j))+dfloat(ie(j))))
rnn=0.5d0*rat*dfloat(ie(j)-is(j))
nn=1+idint(rnn)
is(j)=n_mid-nn
ie(j)=n_mid+nn
15 continue
c
c........................................ Recalculate {wa} according to the
c the new {in1}, {in2} arrays
do 16 i=1,n
wa(i)=1.0d0
16 continue
do 17 j=1,ns
i1=is(j)
i2=ie(j)
do 18 i=i1,i2
wa(i)=0.0d0
18 continue
17 continue
c
return
end
c
c
subroutine solve2_norm_w(m,n1,n2,x,we,xw,c,sig)
c--------------------------------------------------------------------
c This routine solves the Weighted Least Squares problem derived
c from the design matrix {h(i,j)} and data array {x(k)}:
c
c x_est (k) = SUM_j c(j)*h(j,k)
c--------------------------------------------------------------------
c
c DATE: - 2020.02.14
c
c INPUT: - m :: Number of parameters
c - n1 :: First item in {x} where the LS
c solution is to be computed
c - n2 :: Last item in {x} where the LS
c solution is to be computed
c - {xte(j,i)} :: In COMMON block:
c Design matrix
c 1 <= j <= m , 1<= i <= n
c with n1 <= i <= n2
c - {x} :: Data to be fitted
c - {we} :: Data weights
c
c OUTPUT: - {c} :: Regression coefficients
c - {xw} :: Residuals {x-fit}
c - sig :: BIASED sigma
c
c NOTE: - The time indices {i} of the design matrix {h} may
c exceed [n1,n2] in both directions. In the LS
c solution we utilize values only in [n1,n2]
c - SUM{we}=n2-n1+1 MUST be satisfied !
c
c==================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),xw(30333),we(30333),c(995)
dimension g(995,996)
c
common /aa/ xte(995,30333)
c
c------------------------------------------------------------------
c
m1=m+1
c
c............................... Compute normal matrix
do 3 i=1,m
do 4 j=i,m
s=0.0d0
do 5 k=n1,n2
s=s+we(k)*xte(i,k)*xte(j,k)
5 continue
g(i,j)=s
g(j,i)=s
4 continue
3 continue
c
c............................... Compute RHS
do 7 i=1,m
s=0.0d0
do 15 k=n1,n2
s=s+we(k)*xte(i,k)*x(k)
15 continue
g(i,m1)=s
7 continue
c
c............................... Compute solution
call gelim1(g,c,m,IFLAG)
c
c............................... Compute residual {xw} and RMS
sum=0.0d0
do 13 i=n1,n2
s=0.0d0
do 14 j=1,m
s=s+c(j)*xte(j,i)
14 continue
d=x(i)-s
xw(i)=d
sum=sum+we(i)*d*d
13 continue
sig=dsqrt(sum)
c
return
end
c
c
subroutine robust_tfa_four(iwa,n,t,x,m_four,fourf,
&fourts,tfats,we,sig)
c
c*****************************************************************************
c This ROUTINE COMPUTES TFA+FOURIER+AR(edge)-FILTERED TIME SERIES *
c*****************************************************************************
c
c DATE: -- 2020.06.03
c
c INPUT: iwa -- Weighting method key (see MAIN)
c n -- Number of data points of the target time series
c t(i) -- Time values of the target time series
c x(i) -- Values of the target time series
c m_four -- Number of Fourier components
c fourf(i) -- Frequency of the i-th Fourier component
c
c
c OUTPUT: we(i) -- Weights (iterative determined CAUCHY kernel)
c fourts(i) -- Fourier signal computed from the Fourier part
c of the joint TFA-FOUR fit.
c tfats(i) -- Template function derived from the joint
c TFA-FOUR fit
c sig -- *BIASED* estimate of the standard deviation
c of the residuals of the time series with the
c TFA+FOUR components subtracted from the input
c time series.
c
c NOTES: -- nt = ntemp MUST BE satisfied! [ i.e., input time series
c {ta,xa} MUST BE on the TFA TEMPLATE TIMEBASE - see COMMON
c block /tt/ ]
c -- ZERO POINT value is included in the TFA fit
c -- TFA parameters and arrays are pre-computed, and given
c in the COMMON fields
c
c------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 idt(995)
c
dimension t(30333),x(30333),wa(30333),we(30333),xw(30333)
dimension c(995)
dimension fourts(30333),tfats(30333)
dimension fourf(30333),si(30333)
c
common /aa/ xte(995,30333)
common /cc/ m3_tfa,mte0
common /cc1/ idt
common /tt/ ttemp(30333),ntemp
common /rc/ m3_rc
c
c------------------------------------------------------------------------
c
c.................................. Set iteration and other parameters
s3 = 1.0d0
nit = 20
sig0 = 1.0d20
eps = 0.01d0
pi2 = 8.0d0*datan(1.0d0)
c
c s3 = Constant in the weight function {we} is given by s3*sig0,
c where sig0 is the RMS of the time series
c nit = Number of iterations on the weights {we}
c eps = Iteration stops when the relative change in the RMS is
c smaller than eps.
c
c
c.................................. 1 - n0 : TFA templates
c n2 - n3 : FOURIER
m = m_four
m2 = 2*m
n0 = m3_tfa
n2 = n0+1
n3 = n0+m2
c
c.................................. Initialize weights
do 1 k=1,n
wa(k)=1.0d0
we(k)=1.0d0
1 continue
c
c EXTEND TEMPLATE SET BY THE FOURIER COMPONENTS
c 1, sin(om1), cos(om1), sin(om2), cos(om2), ...
c
kk = 0
rn = 1.0d0/dfloat(n)
do 3 k=n2,n3,2
kk = kk+1
omegak = pi2*fourf(kk)
k1 = k
k2 = k+1
ss=0.0d0
sc=0.0d0
do 2 i=1,n
phase = omegak*t(i)
s1=dsin(phase)
c1=dcos(phase)
xte(k1,i)=s1
xte(k2,i)=c1
ss=ss+s1
sc=sc+c1
2 continue
ss=ss*rn
sc=sc*rn
do 4 i=1,n
xte(k1,i)=xte(k1,i)-ss
xte(k2,i)=xte(k2,i)-sc
4 continue
3 continue
m3 = n3
m4 = m3+1
m3_rc= m3
c
c................................ Initialize solution
call solve2_norm_w(m3,1,n,x,we,xw,c,sig)
if(iwa.eq.1) call slide_rms(m_four,n,t,xw,si,wa)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Iteration on {we} STARTS
c
write(*,101)
101 format(/,'>>>>> Robust TFA-FOURIER fit ...')
c
ita=0
do 100 it=1,nit
c
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 100
sig0=sig
ita=ita+1
c................................ Update weights
cc=s3*sig0
c2=cc*cc
s=0.0d0
do 5 i=1,n
dd=xw(i)
d2=dd*dd
d2=1.0d0/(c2+d2)
d2=d2*wa(i)
we(i)=d2
s=s+d2
5 continue
s=1.0d0/s
do 6 i=1,n
we(i)=s*we(i)
6 continue
call solve2_norm_w(m3,1,n,x,we,xw,c,sig)
if(iwa.eq.1) call slide_rms(m_four,n,t,xw,si,wa)
c
100 continue
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Iteration on {we} ENDS
c
c... Compute TFA filter
c
do 26 i=1,n
s=0.0d0
do 27 j=1,n0
s=s+c(j)*xte(j,i)
27 continue
tfats(i)=s
26 continue
c
c... Compute best fitting FOUR time series
c
do 611 i=1,n
s=0.0d0
do 612 j=n2,m3
s=s+c(j)*xte(j,i)
612 continue
fourts(i)=s
611 continue
c
write(*,7) ita
7 format('Actual iteration number in *robust_tfa_four* :',i3)
c
return
end
c
c
subroutine match_obj(sname,n_ref,ref_id,ref_mag,avmag)
c---------------------------------------------------------------------
c This routine matches *sname* with the item in *ref.pos* and
c outputs the corresponding average magnitude *avmag*
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 sname,ref_id(30333)
c
dimension ref_mag(30333)
c
k=0
do 7 i=1,n_ref
if(ref_id(i).eq.sname) k=i
7 continue
if(k.eq.0) write(*,6)
if(k.eq.0) write(*,*) 'Missing reference name =',sname
if(k.eq.0) stop
6 format('No match between *sname* and *ref.pos* !!')
avmag = ref_mag(k)
c
return
end
c
c
subroutine fdec_fast(n0,n,t,x,m,rms)
c---------------------------------------------------------------------
c This yields RMS of {n,t,x} fitted by a FOURIER sum of order *m*
c on the equidistant timebase {t}. The frequencies are multiples
c of f0=1/(t(n)-t(1))
c=====================================================================
c
c n0 = n0 data points are omitted at both edges
c when computing the RMS
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333),y(30333)
dimension s(456),c(456)
c
pi2=8.0d0*datan(1.0d0)
f0=1.0d0/(t(n)-t(1))
dt=(t(n)-t(1))/dfloat(n-1)
rn2=2.0d0/dfloat(n)
c
i1=n0+1
i2=n-n0
c
c0=0.0d0
do 1 i=1,n
c0=c0+x(i)
1 continue
c0=c0/dfloat(n)
do 8 i=1,n
y(i)=x(i)-c0
8 continue
if(m.gt.0) go to 4
c
c.................................. Compute RMS for m = 0
disp=0.0d0
do 5 i=i1,i2
disp=disp+y(i)*y(i)
5 continue
go to 6
4 continue
c
c.................................. Compute RMS for m > 0
do 2 j=1,m
wj=pi2*f0*dfloat(j)*dt
s1=0.0d0
c1=0.0d0
do 3 i=1,n
alpha=wj*dfloat(i-1)
xx=y(i)
s1=s1+xx*dsin(alpha)
c1=c1+xx*dcos(alpha)
3 continue
s(j)=s1*rn2
c(j)=c1*rn2
2 continue
c
disp=0.0d0
w0=pi2*f0*dt
do 2000 i=i1,i2
sum=0.0d0
w=w0*dfloat(i-1)
do 2001 j=1,m
alpha=w*dfloat(j)
sum=sum+s(j)*dsin(alpha)+c(j)*dcos(alpha)
2001 continue
devi=y(i)-sum
disp=disp+devi*devi
2000 continue
c
6 continue
rms=dsqrt(disp/dfloat(i2-i1+1-2*m-1))
c
return
end
c
c
subroutine four_order_fast(isurvey,ext,n,t,x,m_ord)
c---------------------------------------------------------------------
c This is to find the OPTIMUM FOURIER order to fit {n,t,x}
c---------------------------------------------------------------------
c
c DATE: 2020.06.03
c
c INPUT : {n,t,x} = Time series to be tested
c ext = In computing RMS, ext*n data points at the
c edges are omited
c isurvey = If 1, than no output file "four_ord.dat"
c
c OUTPUT: m_ord = Optimum order of the Fourier series
c
c---------------------------------------------------------------------
c
c NOTES: - Data are transformed to equidistant timebase, then
c simple DFT is run on this dataset, yielding the exact
c LS solution (due to the equidistant timebase).
c - Frequency sets tested: f0, (f0,2*f0), (f0,2*f0,3*f0),
c ... (f0, 2*f0, ..., 90*f0), with f0=1/(t(n)-t(1))
c - It is assumed that the time series is ordered, i.e.,
c t(i+1) > t(i)
c - The optimum Fourier order is determined in the following
c way:
c 1. Compute UNBIASED sigma_fit values for a sequence
c of Fourier orders from m1 to m2 with step ms
c 2. Robustly fit a straight line to the last 2/3 -rd
c of the sigma_fit values
c 3. Compute the UNWEIGHTED UNBIASED RMS of the fit
c (parameter *stdv*)
c 4. Extrapolate this line all the way to the lowest
c Fourier order
c 5. Compute the difference between the actual RMS of the
c Fourier fit and the extrapolated value from the
c lowest Fourier order to the highest one
c 6. The 1st kind of optimum Fourier order is defined at
c the lowest order at which the difference is still
c positive between the measured and extrapolted value,
c and it is larger than *stdv*
c 7. The 2nd kind of optimum Fourier order is derived
c on the basis of the relative decrease of the variances
c at different Fourier order.
c 8. The FINAL OPTIMUM value of the Fourier order comes
c from the combination of the 1st and 2nd kind optimum
c values and a minimum order condition.
c [see the paper for more details]
c
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),u(30333),v(30333),x(30333),xi(30333)
dimension xd(30333),yd(30333),z(30333),r2(30333)
dimension si(456)
dimension mm(456)
c
n0 = idint(ext*dfloat(n))
m1 = 1
m2 = 150
ms = 3
f23= 2.0d0/3.0d0
r2m= 6.0d0
c
c.................................... Compute average sampling rate
t1=t(1)
dt=(t(n)-t1)/dfloat(n-1)
c
c.................................... Interpolate data on that grid
IQ=1
N1=n
N2=n
do 6 i=1,n
u(i)=t(i)
v(i)=x(i)
xd(i)=t1+dt*dfloat(i-1)
6 continue
CALL INTER2(IQ,u,v,N1,xd,yd,N2)
c
c.................................... START cycle of testing the FOURIER ORDER
k=1
m=0
call fdec_fast(n0,n,xd,yd,m,sig)
mm(k)=m
si(k)=sig
c
do 2 m=m1,m2,ms
k=k+1
call fdec_fast(n0,n,xd,yd,m,sig)
si(k)=sig
mm(k)=m
2 continue
c
c... Compute sigma(var) for the last Fourier order
c
s_va2=dsqrt(2.0d0)*si(k)*si(k)/dsqrt(dfloat(n-2*n0))
c
c... Robust linear fit to the high-m fraction part of the
c m --> sigma function
c
i1=idint(f23*dfloat(k))
i2=k
nn=0
np=1
do 7 i=i1,i2
nn=nn+1
u(nn)=dfloat(mm(i))
v(nn)=si(i)
7 continue
call robust_pol(nn,np,u,v,xi,sig,stdv)
c
if(isurvey.eq.0) open(28,file='four_ord.dat')
c
c............................. Extrapolation of the linear fit
do 9 j=1,nn
i=k-j+1
j1=nn-j+1
u(i)=u(j1)
v(i)=v(j1)
xi(i)=xi(j1)
9 continue
u1=u(i1)
u2=u(i2)
v1=xi(i1)
v2=xi(i2)
rr=(v2-v1)/(u2-u1)
do 8 i=1,k
if(i.lt.i1) z(i)=v1+(mm(i)-u1)*rr
if(i.ge.i1) z(i)=xi(i)
rat1=(si(i)-z(i))/sig
r2(i)=si(i)*si(i)/s_va2
if(isurvey.eq.0) write(28,10) mm(i),si(i),z(i),rat1,r2(i)
10 format(i4,2(1x,f10.8),2(1x,f10.1))
8 continue
close(28)
c
c... Find the last mm(i) with si(i)-z(i) > sdev
c
k0=0
k1=0
rr=r2(k)
do 4 i=1,k
if(k0.gt.0) go to 5
if(si(i)-z(i).gt.stdv) go to 5
k0=i
5 continue
if(k1.gt.0) go to 4
if(r2(i)-rr.gt.r2m) go to 4
k1=i
4 continue
m_ord = mm(k0)
if(mm(k1).lt.m_ord) m_ord=mm(k1)
m_ord=m_ord+10
c
write(*,*) 'k0, k1, mm(k0), mm(k1) =',k0,k1,mm(k0),mm(k1)
c
if(m_ord.lt.20) m_ord=20
c
return
end
c
c
subroutine ar_vs_four(npr,n,z,xf,xfa,s_xf1,s_xfa1,
&s_xf2,s_xfa2,iflag1,iflag2)
c-------------------------------------------------------------------------
c This routine compares the clipped RMS values for the AR(xfa) and
c FOUR(fourts) approximations at the edges, and modifies {xfa} if
c FOUR proves to show lower RMS.
c=========================================================================
c
c INPUT: npr - Datapoint number at the edges for the AR modeling
c n - Total number of datapoints
c {z} - Systematics-free time series: {z}={x_int}-{tfats}
c {xf} - Fourier fit
c {xfa} - Fourier fit with AR/polynomial edge approximation
c
c OUTPUT: s_xf1 - RMS of the residuals by FOUR at the LEFT edge
c s_xfa1- RMS of the residuals by AR at the LEFT edge
c s_xf2 - RMS of the residuals by FOUR at the RIGHT edge
c s_xfa2- RMS of the residuals by AR at the RIGHT edge
c iflag - 1 - FOUR yields lower residual
c 2 - AR yields lower residual
c {xfa} - Updated {xfa}
c
c-------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension xfa(30333),z(30333),xf(30333)
dimension xf1(30333),xf2(30333)
dimension xfa1(30333),xfa2(30333)
c
arsig=3.0d0
c
k1=0
k2=0
do 5 i=1,n
xx=z(i)
if(i.gt.npr) go to 6
k1=k1+1
xf1(k1)=xx-xf(i)
xfa1(k1)=xx-xfa(i)
go to 5
6 continue
if(i.lt.n-npr) go to 5
k2=k2+1
xf2(k2)=xx-xf(i)
xfa2(k2)=xx-xfa(i)
5 continue
c
call robust_aver(k1,xf1,a_xf1,s_xf1,sigd)
call robust_aver(k1,xfa1,a_xfa1,s_xfa1,sigd)
call robust_aver(k2,xf2,a_xf2,s_xf2,sigd)
call robust_aver(k2,xfa2,a_xfa2,s_xfa2,sigd)
c
iflag1=1
if(s_xfa1.lt.s_xf1) iflag1=2
iflag2=1
if(s_xfa2.lt.s_xf2) iflag2=2
if(iflag1.eq.2) go to 15
do 4 i=1,npr
xfa(i)=xf(i)
4 continue
15 continue
if(iflag2.eq.2) go to 16
do 3 i=1,npr
j=n-npr+i
xfa(j)=xf(j)
3 continue
16 continue
c
write(*,17) s_xf1,s_xfa1,iflag1
17 format('Weighted s_xf1,s_xfa1,iflag1 = ',2(1x,f8.6),1x,i2)
write(*,18) s_xf2,s_xfa2,iflag2
18 format('Weighted s_xf2,s_xfa2,iflag2 = ',2(1x,f8.6),1x,i2)
if(iflag1.eq.1) write(*,*) 'LEFT: FOUR'
if(iflag1.eq.2) write(*,*) 'LEFT: AR'
if(iflag2.eq.1) write(*,*) 'RIGHT: FOUR'
if(iflag2.eq.2) write(*,*) 'RIGHT: AR'
c
return
end
c
c
subroutine robust_aver(n,x,ave,sig,sigd)
c-----------------------------------------------------------------
c Computes robust average and RMS
c=================================================================
c
c NOTE: The constant in the weight factor is given by:
c (s3*RMS)**2, where RMS is the RMS of the past fit
c-----------------------------------------------------------------
c
c INPUT: n = Number of data points
c {x} = Target array to process
c
c OUTPUT: ave = Robust average
c sig = Weighted standard deviation of {x}
c sigd = Equal-weighted standard deviation of {x}
c
c==============================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension x(30333),we(30333)
c
c================================== For TESTING
c n=400
c iseed=2000*136+1
c si=0.83d0
c do 11 i=1,n
c x(i)=2.0d0+si*gauss(iseed)
c 11 continue
c do 12 i=10,40
c x(i)=0.0d0+si*gauss(iseed)
c x(n-i)=-1.0d0+si*gauss(iseed)
c 12 continue
c
c do 13 i=1,n
c write(33,*) i,x(i)
c 13 continue
c================================== For TESTING
c
c avsig= 2.5d0
c sig3 = 4.0d0
c
s3 = 1.0d0
nit = 20
sig = 1.0d30
sig0 = 1.0d0
eps = 0.001d0
call astat(x,n,ave,sig0)
c
c D = SUM (x-f)^2/(cc+(x-f)^2)
c
c................................ START ITERATION CYCLE
do 1 it=1,nit
c
if(it.le.2) go to 4
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 1
sig0=sig
4 continue
c
c................................ Update weights
cc=s3*sig0
c2=cc*cc
s=0.0d0
do 5 i=1,n
dd=x(i)-ave
d2=dd*dd
d2=1.0d0/(c2+d2)
s=s+d2
we(i)=d2
5 continue
s=1.0d0/s
do 6 i=1,n
we(i)=s*we(i)
6 continue
c
s=0.0d0
do 9 i=1,n
s=s+we(i)*x(i)
9 continue
ave=s
s=0.0d0
sd=0.0d0
do 8 i=1,n
d=x(i)-ave
d2=d*d
sd=sd+d2
s=s+we(i)*d2
8 continue
sig=dsqrt(s)
sigd=dsqrt(sd/dfloat(n-1))
1 continue
c
return
end
c
c
subroutine single_tran(n,t,x,tmin,tmax,dsp_min,
&tce_1,dep_1,qtr_1,qin_1,x_1,nev_1,ndit_1,itoc_1,dsp_1,nsing)
c---------------------------------------------------------------------
c Search for single transit event
c---------------------------------------------------------------------
c
c INPUT: - dsp_min= Minimum DSP for the SINGLE transit
c qualification
c - n = Number of data points
c - {t} = Moments of time
c - {x} = Time series values
c - tmin = MIN({t})
c - tmax = MAX({t})
c
c OUTPUT: - tce_1 = Moment of the center of the transit
c - dep_1 = Depth of the transit
c - qtr_1 = t14/P_orb
c - qin_1 = t12/t14
c - {x_1} = Best fit to the single transit
c - nev_1 = Number of transit events
c - ndit_1 = Number of data points in the transit
c - itoc_1 = Transit occupation code
c - dsp_1 = Dip Significance Parameter
c - nsing = 0 - The transit found is NOT a single one
c = 1 - The transit found IS a SINGLE one, based
c on nev_1 and DSP
c
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension u(30333),v(30333)
dimension xft(30333),x_1(30333)
dimension p(155555)
dimension intr(30333)
c
nb = 1000
qmi = 0.002
qma = 0.030
nf = 1000
tot = tmax-tmin
fmin = 1.0d0/tot
totdf = 2.0d0/tot
fmax = fmin+totdf
df = (fmax-fmin)/dfloat(nf-1)
c
call eebls(n,t,x,u,v,nf,fmin,df,nb,qmi,qma,
&p,bper,bpow,depth,qtran,in1,in2)
c
p0=bper
call foldts4_new(p0,qmi,qma,nb,n,t,x,xft,
&intr,tdepth,ring,rtran,tcen,nintr,sig4)
c
f0=1.0d0/bper
call events(n,t,f0,tcen,qtran,nev,ndit,itoc)
phas=(tcen-tmin)/tot
c
call dsp_ok(n,t,x,p0,tcen,qtran,dsp,snrt,fret)
c
tce_1 = tcen
dep_1 = tdepth
qtr_1 = rtran
qin_1 = ring
do 2 i=1,n
x_1(i) = xft(i)
2 continue
c
c... Checking SINGLE transit criteria
c
nsing = 0
if(nev.eq.1.and.dsp.gt.dsp_min) nsing=1
c
write(*,1) p0/tot,phas,tdepth,nev,ndit,itoc,dsp
1 format(/,38('-'),/,
&'SINGLE transit test:',/,
&'Period/total time span = ',f8.3,/,
&'Phase of T_cen = ',f8.3,/,
&'Transit depth = ',f8.6,/,
&'Number of transit events = ',i8,/,
&'Number of intransit points = ',i8,/,
&'Transit occupation code = ',i8,/,
&'DSP = ',f8.3,/,
&38('-'))
c
tce_1 = tcen
dep_1 = tdepth
qtr_1 = rtran
qin_1 = ring
nev_1 = nev
ndit_1= ndit
itoc_1= itoc
dsp_1 = dsp
c
if(nsing.eq.0) go to 3
write(*,4)
4 format(/,40('*'),/,
&'* This target contains a SINGLE TRANSIT *',/,
&40('*'),/)
3 continue
c
return
end
c
c
subroutine multi_rec(ksm,nt,ta,xi,xa,nbb,
&qmi1,qma1,ntw,fw,dw,snrw,snr_min,asig2,
&ntc,xtr,xft,ft,xt,we,fr,td,qt,ri,tc,sig,nclip)
c
c
c >>>>> The OUTPUT of the ksm=0 case should still be fixed !!!!
c
c
c
c
c***************************************************************************
c MULTI-PERIODIC TRANSIT SIGNAL RECONSTRUCTION *
c***************************************************************************
c
c INPUT: ksm -- 0 = no reconstruction; 1 = reconstruction
c nt -- Number of data points
c {ta} -- Timebase ( {ta}={ttemp} )
c {xi} -- Original time series, interpolated to the TFA
c template timebase ( {xi}={x_int} )
c {xa} -- Some systematics-filtered version of the target
c time series [to compute initial estimate on the
c transit signal {xft}; {xa}={x}]
c {we} -- Weights from the earlier dominant transit
c estimates (via *one_rec*)
c nbb -- Number of bins for "eebls" in "foldts4_w"
c qmi1 -- Minimum q_tran for "eebls" in "foldts4_w"
c qma1 -- Maximum q_tran for "eebls" in "foldts4_w"
c ntw -- Number of candidate transit signal components
c {fw} -- Frequencies of the above transit signals
c {dw} -- Transit depths of the above transit signals
c {snrw} -- BLS SNR values of the transit components
c snr_min -- Reconstruction is performed for transit signal
c candidates that satisfy these conditions:
c dw(i) > 0; snrw(i) > snr_min
c
c OUTPUT: ntc -- Number of transit components satisfying the
c dw(i) > 0; snrw(i) > snr_min conditions
c {xtr} -- Outlier filtered (sigma-clipped), systematics-
c and stellar variability- free time series
c {xft} -- Robust composite (i.e., containing all
c significant transit components) fit to {xtr}
c {ft(j,i)} -- Fitted transit signal of the j-th component
c {xt(j,i)} -- Outlier filtered (sigma-clipped), systematics-
c and stellar variability- free transit signal
c of the j-th component
c {we} -- Weights from the best fitting full model
c {fr} -- Frequencies of the output components
c {td} -- Transit depths
c {qt} -- Relative transit durations
c {ri} -- t12/t14
c {tc} -- Transit centers
c sig -- **BIASED** estimate of the standard deviation
c of the residuals of the time series with the
c TFA+FOUR+TRANSIT components subtracted from the
c input time series. Outliers are included in the
c residuals. All items are taken with equal weights.
c nclip -- Number of sigma-clipped outliers in the composite
c time series {xtr}. [Outliers are selected in respect
c of {xft}.]
c
c
c NOTES: -- {xi} = {x_int} in MAIN
c -- {xa} = {x} in MAIN ({xdat} in "datain_tfa_four")
c [i.e., {xa} is the "non-reconstructed" time
c series]
c -- TFA parameters and arrays are pre-computed, and given
c in the COMMON fields
c -- The full reconstruction includes:
c transit + systematics + stellar variability;
c {xte} includes systematics + stellar variability
c -- There are altogether "m3_rc" correcting vectors.
c We add the updated approximation of the transit signal
c {xft} at each iteration to the above set of correcting
c vectors. Thereby we enable the full reconstruction.
c -- This routine uses the transit shapes derived from the
c *non-reconstructed* input signal {xa} as computed in
c the first part of this routine. During the reconstruc-
c tion we perform a WEIGHTED LS fit ONLY to the TRANSIT
c DEPTHS and the shapes are left UNCHANGED, as derived
c in the first part of this routine.
c -- The full model is derived by a weighted Least Squares
c fit to handle outliers.
c -- Output time series are the ones obtained by the above
c fit, but filtered (sigma-clipped) from outliers.
c -- If short-period oscillations prevail, then some of
c the components may go wild, following the pattern of
c the oscillations.
c
c------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension xi(30333),ta(30333),xa(30333),we(30333)
dimension xft(30333),xtr(30333),xw(30333),x(30333),xf(30333)
dimension fw(100),dw(100),snrw(100)
dimension td(100),ri(100),qt(100),tc(100),tp(100),fr(100)
dimension ft(6,30333),xt(6,30333)
dimension g(995,996),c(995)
dimension intr(30333)
c
common /aa/ xte(995,30333)
common /rc/ m3_rc
c
c------------------------------------------------------------------------
c
c... Set iteration parameters, etc.
c
s3 = 1.0d0
nit = 15
ita = 0
nit5 = 5
eps = 0.001d0
sig0 = 1.0d33
m3 = m3_rc+1
m4 = m3+1
m5 = m3_rc
c
write(*,*) '>>>>> ITERATION in *multi_rec*'
c
c#####################################################################
c
c... FIRST estimate of the components [equal weights]
c
k=0
do 201 jt=1,ntw
if(snrw(jt).lt.snr_min) go to 201
if(dw(jt).lt.0.0d0) go to 201
p0=1.0d0/fw(jt)
k=k+1
call foldts4_new(p0,qmi1,qma1,nbb,nt,ta,xa,xft,
&intr,tdepth,ring,rtran,epc,nintr,sig4)
fr(k)=fw(jt)
td(k)=tdepth
ri(k)=ring
qt(k)=rtran
tc(k)=epc
do 200 i=1,nt
ft(k,i)=xft(i)
200 continue
201 continue
ntc=k
c
c... ITERATIVE estimate of the components [equal weights]
c
do 203 it=1,nit5
do 204 k=1,ntc
do 207 j=1,ntc
tp(j)=1.0d0
207 continue
tp(k)=0.0d0
do 205 i=1,nt
s=xa(i)
do 206 j=1,ntc
s=s-tp(j)*ft(j,i)
206 continue
x(i)=s
205 continue
p0=1.0d0/fw(k)
call foldts4_new(p0,qmi1,qma1,nbb,nt,ta,x,xft,
&intr,tdepth,ring,rtran,epc,nintr,sig4)
td(k)=tdepth
ri(k)=ring
qt(k)=rtran
tc(k)=epc
do 208 i=1,nt
ft(k,i)=xft(i)
208 continue
sig=sig4
204 continue
203 continue
c
if(ksm.eq.0) go to 300
c
c#####################################################################
c
m3 = m3_rc+ntc
m4 = m3+1
m5 = m3_rc
k1 = m3_rc+1
k2 = m3
c
c#####################################################################
c... START ITERATION on the weights of the FULL model
c
do 100 it=1,nit
c
rsig=dabs(1.0d0-sig/sig0)
if(rsig.lt.eps) go to 100
ita=ita+1
sig0=sig
c
c... Add transit signals to the TFA template set
c
do 212 j=k1,k2
i=j-k1+1
do 32 k=1,nt
xte(j,k)=ft(i,k)
32 continue
212 continue
c
c... Compute TFA normal matrix
c
do 601 i=1,m3
do 602 j=i,m3
s=0.0d0
do 615 k=1,nt
s=s+we(k)*xte(i,k)*xte(j,k)
615 continue
g(i,j)=s
g(j,i)=s
602 continue
601 continue
c
c... Compute RHS
c
do 303 i=1,m3
s=0.0d0
do 306 k=1,nt
s=s+we(k)*xi(k)*xte(i,k)
306 continue
g(i,m4)=s
303 continue
c
c... Compute LS solution
c
call gelim1(g,c,m3,IFLAG)
c
c... Compute transit signal
c {xtr} = {xi}-{TFA}-{FOUR}
c and standard deviation of the residual
c {res} = {xi}-{TFA}-{FOUR}-{TRANSIT}
c
d=0.0d0
do 26 k=1,nt
s1=0.0d0
do 27 j=1,m5
s1=s1+c(j)*xte(j,k)
27 continue
xx=xi(k)-s1
s2=0.0d0
do 213 j=k1,k2
i=j-k1+1
ftik=c(j)*xte(j,k)
s2=s2+ftik
ft(i,k)=ftik
213 continue
yy=xx-s2
d=d+yy*yy
xtr(k)=xx
xft(k)=s2
xw(k)=yy
26 continue
sig=dsqrt(d/dfloat(nt))
c
c.................................. UPDATE {we}
cc=s3*sig
c2=cc*cc
s=0.0d0
do 5 k=1,nt
dd=xw(k)
d2=dd*dd
d2=1.0d0/(c2+d2)
we(k)=d2
s=s+d2
5 continue
s=1.0d0/s
do 6 k=1,nt
we(k)=s*we(k)
6 continue
c
c-------------------------------------------------------------
c Unfortunately, the transit shape matters in computing
c the transit depth (e.g., 211817229 from C05). So, we
c need to update the individual transits - one-by-one ...
c-------------------------------------------------------------
c
c.................................. UPDATE {xt}, {ft}
do 40 i=1,nt
x(i)=xtr(i)-xft(i)
40 continue
do 41 j=1,ntc
do 42 i=1,nt
xt(j,i)=x(i)+ft(j,i)
42 continue
41 continue
c
do 43 j=1,ntc
do 44 i=1,nt
x(i)=xt(j,i)
44 continue
p0=1.0d0/fr(j)
call foldts4_w(p0,qmi1,qma1,nbb,nt,ta,x,xf,
&we,intr,tdepth,ring,rtran,epc,nintr,sig4)
do 45 i=1,nt
ft(j,i)=xf(i)
45 continue
43 continue
c
100 continue
c
write(*,1) ita
1 format('Actual iteration number in *multi_rec* :',i3)
c
c
c... END ITERATION on {ft} for the FULL model
c#####################################################################
c
300 continue
c
c... Prepare output arrays
c ntc, {xtr}, {xft}, {xt}, {ft}, {we}, {td}, {qt}, {ri}, {tc}
c
do 223 i=1,nt
x(i)=xtr(i)-xft(i)
223 continue
do 222 j=1,ntc
do 224 i=1,nt
xt(j,i)=x(i)+ft(j,i)
224 continue
222 continue
c
c... Sigma-clipping of the composite time series {xtr}, {xt}
c
call sigma_clipping_r(-1,asig2,nt,xtr,xft,nclip)
c
c... Sigma-clipping of the individual transit time series {xt}
c
do 216 j=1,ntc
do 217 i=1,nt
x(i)=xt(j,i)
xf(i)=ft(j,i)
217 continue
call sigma_clipping_r(-1,asig2,nt,x,xf,ncl)
do 219 i=1,nt
xt(j,i)=x(i)
219 continue
216 continue
c
c... Compute transit parameters and save densely sampled transit fits
c
do 220 j=1,ntc
do 221 i=1,nt
x(i)=ft(j,i)
221 continue
p0=1.0d0/fr(j)
call foldts4_new(p0,qmi1,qma1,nbb,nt,ta,x,xf,
&intr,tdepth,ring,rtran,tcen,nintr,sig4)
td(j)=tdepth
qt(j)=rtran
ri(j)=ring
tc(j)=tcen
ktr=j
d00=tdepth
q00=rtran
r00=ring
call save_trc(ktr,d00,q00,r00)
220 continue
c
return
end
c
c
subroutine save_trc(ktr,d00,q00,r00)
c--------------------------------------------------------------------------
c This routine saves the densely sampled neighborhood of a transit fit
c--------------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*10 fname
dimension u(30333),tru(30333)
c
c.................................. Set parameters for the densely-sampled fit
n_ph = 200
r_q = 1.5d0
c
c.................................. Set output file name
if(ktr.eq.1) fname='tr1_fit.lc'
if(ktr.eq.2) fname='tr2_fit.lc'
if(ktr.eq.3) fname='tr3_fit.lc'
if(ktr.eq.4) fname='tr4_fit.lc'
if(ktr.eq.5) fname='tr5_fit.lc'
if(ktr.eq.6) fname='tr6_fit.lc'
c
c... Generate densely sampled transit fit
c
p00=1.0d0
e00=0.0d0
ph_min=-r_q*q00
ph_max= r_q*q00
d_ph=(ph_max-ph_min)/dfloat(n_ph-1)
do 1 i=1,n_ph
u(i)=ph_min+d_ph*dfloat(i-1)
1 continue
call trapu(n_ph,u,p00,e00,q00,r00,tru)
open(7,file=fname)
u00=-0.5d0
v00=0.0d0
write(7,3) u00,v00
do 2 i=1,n_ph
write(7,3) u(i),d00*tru(i)
3 format(f7.4,2x,f9.6)
2 continue
u00=0.5d0
write(7,3) u00,v00
close(7)
c
return
end
c
c
subroutine save_m_rec(sname,ntc,n,t,xtr,xft,ft,xt,
&fr,td,qt,ri,tc)
c---------------------------------------------------------------------
c This routine saves the reconstructed time series (un-folded and
c folded for the individual components) and the transit parameters
c---------------------------------------------------------------------
c NOTES: - We set the OOT level to zero for each transit component
c=====================================================================
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*6 f(6)
character*9 sname
c
dimension t(30333),xft(30333),xtr(30333),x(30333)
dimension ft(6,30333),xt(6,30333)
dimension td(100),ri(100),qt(100),tc(100),fr(100)
c
c... This is a primitive setting the folded LC file names
c
f(1)='tr1.lc'
f(2)='tr2.lc'
f(3)='tr3.lc'
f(4)='tr4.lc'
f(5)='tr5.lc'
f(6)='tr6.lc'
c
c... Save phase-folded LCs for the individual transit components
c
c0=0.0d0
c
do 2 j=1,ntc
open(1,file=f(j))
f0=fr(j)
epoch=tc(j)
do 1 i=1,n
ph = f0*(t(i)-epoch)
ph = ph - idint(ph)
if(ph.gt. 0.5d0) ph=ph-1.0d0
if(ph.lt.-0.5d0) ph=ph+1.0d0
write(1,3) ph,xt(j,i),ft(j,i)
3 format(3(1x,f9.6))
1 continue
close(1)
2 continue
c
c... Save un-folded LC and composite fit containing all transits
c
open(1,file='tr_all.lc')
do 7 i=1,n
write(1,8) t(i),xtr(i),xft(i)
8 format(f13.6,2(2x,f9.6))
7 continue
close(1)
c
c... Save transit parameters
c
open(25,file='multi_par.dat')
write(25,38)
write(*,38)
38 format('#=====================================================',
&'===================',/,
&'# EPIC f [c/d] depth t14/P t12/t14 T_c - ',
&'2400000 DSP',/,
&'#------------------------------------------------------------',
&'------------')
c
do 4 j=1,ntc
do 6 i=1,n
x(i)=xt(j,i)
6 continue
p0=1.0d0/fr(j)
tcen=tc(j)
depth=td(j)
qtran=qt(j)
call dsp_ok(n,t,x,p0,tcen,qtran,dsp,snrt,fret)
if(dsp.gt.99.9d0) dsp=99.9d0
write(25,5) sname,fr(j),td(j),qt(j),ri(j),tc(j),dsp
write(*,5) sname,fr(j),td(j),qt(j),ri(j),tc(j),dsp
5 format(a9,2x,f11.7,2x,f8.6,2(2x,f7.5),2x,f15.6,1x,f5.1)
4 continue
close(25)
c
return
end
C
C
subroutine gelim1(A,X,M,IFLAG)
C
C SUBROUTINE FOR SOLVING SIMULTANEOUS LINEAR EQUATIONS BY GAUSSIAN
C ELIMINATION WITH PARTIAL PIVOTING.
C
c=================================================================
c
c A Coefficient matrix A extended to include vector B as its
c last column
c X Solution vector X
c M Size of the matrix A
c IFLAG Integer flag for a singular or very ill-conditioned
c coefficient matrix (normally returned as 0 to the calling
c program, and 1 if the ill-conditioning test is satisfied)
c
c=================================================================
c From: https://www.crcpress.com/downloads/K12712/WEB_Material.pdf
c Date: 2020.03.27
c Note: Slightly simplified, because of the standard format of the
c normal equations used to solve the given LS problem
c-----------------------------------------------------------------
c
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER*4 (I-N)
C
DIMENSION A(995,996),X(995)
C
EPS=1.0D-50
NEQN=M
C
C INITIALIZE ILL-CONDITIONING FLAG.
IFLAG=0
C
C SCALE EACH EQUATION TO HAVE A MAXIMUM COEFFICIENT MAGNITUDE OF UNITY.
JMAX=NEQN+1
DO 2 I=1,NEQN
AMAX=0.0d0
DO 1 J=1,NEQN
AMAX=DMAX1(AMAX,DABS(A(I,J)))
1 CONTINUE
RAMAX=1.0D0/AMAX
DO 12 J=1,JMAX
A(I,J)=A(I,J)*RAMAX
12 CONTINUE
2 CONTINUE
C
C COMMENCE ELIMINATION PROCESS.
DO 5 K=1,NEQN-1
C
C SEARCH LEADING COLUMN OF THE COEFFICIENT MATRIX FROM THE DIAGONAL
C DOWNWARDS FOR THE LARGEST VALUE.
IMAX=K
DO 3 I=K+1,NEQN
IF(DABS(A(I,K)).GT.DABS(A(IMAX,K))) IMAX=I
3 CONTINUE
C
C IF NECESSARY, INTERCHANGE EQUATIONS TO MAKE THE LARGEST COEFFICIENT
C BECOME THE PIVOTAL COEFFICIENT.
IF(IMAX.NE.K) THEN
DO 4 J=K,JMAX
ATEMP=A(K,J)
A(K,J)=A(IMAX,J)
A(IMAX,J)=ATEMP
4 CONTINUE
END IF
C ELIMINATE X(K) FROM EQUATIONS (K+1) TO NEQN, FIRST TESTING FOR
C EXCESSIVELY SMALL PIVOTAL COEFFICIENT (ASSOCIATED WITH A SINGULAR
C OR VERY ILL-CONDITIONED MATRIX).
IF(DABS(A(K,K)).LT.EPS) THEN
IFLAG=1
RETURN
END IF
RAKK=1.0d0/A(K,K)
DO 25 I=K+1,NEQN
FACT=A(I,K)*RAKK
DO 15 J=K,JMAX
A(I,J)=A(I,J)-FACT*A(K,J)
15 CONTINUE
25 CONTINUE
5 CONTINUE
C SOLVE THE EQUATIONS BY BACK SUBSTITUTION, FIRST TESTING
C FOR AN EXCESSIVELY SMALL LAST DIAGONAL COEFFICIENT.
IF(DABS(A(NEQN,NEQN)).LT.EPS) THEN
IFLAG=1
RETURN
END IF
X(NEQN)=A(NEQN,JMAX)/A(NEQN,NEQN)
DO 7 I=NEQN-1,1,-1
SUM=A(I,JMAX)
DO 6 J=I+1,NEQN
SUM=SUM-A(I,J)*X(J)
6 CONTINUE
X(I)=SUM/A(I,I)
7 CONTINUE
C
RETURN
END
C
C
subroutine check_sys_four(n,t,x,f03,snr03,f35,snr35)
c----------------------------------------------------------------
c This routine performs fast DFT on the systematic- stellar
c variability-free data to check if there are still traces of
c these unwanted parts of the original signal
c----------------------------------------------------------------
c INPUT: - {t(i),x(i); i=1,2,...,n} Input time series
c
c OUTPUT: - f03 :: Peak frequency in [0,3] c/d
c - snr03 :: SNR of f03
c - f35 :: Peak frequency in [3,5] c/d
c - snr35 :: SNR of f35
c
c NOTES: - The [0,3] frequency interval is relevant for those
c variables that are not captured by the blind
c Fourier fit employed in this code (thought to be
c effective in ~[0,1]c/d.
c - For the systematics we focus in the [3,5] c/d
c frequency interval centered around the ~4c/d
c thruster ignation cycle, characteristics of the K2
c mission.
c================================================================
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension t(30333),x(30333)
dimension p(155555)
c
nf=2500
f1=0.0d0
f2=3.0d0
df=(f2-f1)/dfloat(nf-1)
call fdft(n,t,x,nf,f1,df,p,pfr,pam)
call spsnr1(nf,p,f1,df,f03,a03,snr03,spd03)
f3=3.0d0
f4=5.0d0
nf=1500
df=(f4-f3)/dfloat(nf-1)
call fdft(n,t,x,nf,f3,df,p,pfr,pam)
call spsnr1(nf,p,f3,df,f35,a35,snr35,spd35)
c
write(*,1) f1,f2,f03,a03,snr03,f3,f4,f35,a35,snr35
1 format(/,'Main DFT peak in the TFA+FOUR-filtered data:',/,
&' fmin fmax f_peak a_peak SNR_peak',/,
&2(1x,f4.1),2(1x,f10.6),1x,f5.1,/,
&2(1x,f4.1),2(1x,f10.6),1x,f5.1)
c
return
end
c
c
subroutine four_vs_tfa(n,tfats,xfa,sig_four,sig_tfa)
c---------------------------------------------------------------------
c This routine computes RMS of {tfats} and {xfa}
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
dimension tfats(30333),xfa(30333)
c
call astat(tfats,n,aver,sig_tfa)
call astat(xfa,n,aver,sig_four)
c
return
end
c
c
subroutine save_pre_bls(isflc,sname,n,t,x,xfa,tfats,we,x_rec)
c---------------------------------------------------------------------
c This routine saves pre-BLS time series
c---------------------------------------------------------------------
c
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
c
character*9 sname
character*16 fn_save
dimension t(30333),x(30333),xfa(30333),tfats(30333),we(30333)
dimension x_rec(30333)
c
c....................................... Save pre-bls LC
open(27,file='pre_bls.dat')
do 2 i=1,n
x_rec(i) = x(i)
write(27,3) t(i),x(i),xfa(i),we(i)
3 format(f13.5,3(2x,f11.6))
2 continue
close(27)
c
if(isflc.eq.0) go to 1
c
c........................................ Save LC as flc/sname.lc
fn_save='flc/'//sname//'.lc'
open(31,file=fn_save)
do 4 i=1,n
write(31,3) t(i),x(i),xfa(i),tfats(i)
4 continue
close(31)
c
1 continue
c
return
end
c
c