| 1 |
PRO DOWEED, infile, outfile, icasea=icasea, ndt=ndt, tsample=tsample, weedthresh=weedthresh, prev=prev |
| 2 |
|
| 3 |
if (n_elements(ndt) eq 0) then ndt=103680 |
| 4 |
if (n_elements(tsample) eq 0) then tsample=60.0 |
| 5 |
;the above works for any 72d timeseries, since only the product of ndt and tsample is used |
| 6 |
if (n_elements(icasea) eq 0) then icasea=2 |
| 7 |
if (n_elements(weedthresh) eq 0) then weedthresh=10 |
| 8 |
|
| 9 |
; get rid of following, now pass all arguments |
| 10 |
; parms.idl still exists, used by scripts |
| 11 |
; parms.idl should define ndt, tsample, and weedthresh |
| 12 |
;spawn,'cat parms.idl',hold |
| 13 |
;for i=0,n_elements(hold)-1 do idlisstupid=execute(hold[i]) |
| 14 |
|
| 15 |
binmics=1e6/(ndt*tsample) |
| 16 |
|
| 17 |
;numult=8 |
| 18 |
;ampmult=8 |
| 19 |
;widmult=8 |
| 20 |
;amult=8 |
| 21 |
;asmult=8 |
| 22 |
|
| 23 |
numult=weedthresh |
| 24 |
ampmult=weedthresh |
| 25 |
widmult=weedthresh |
| 26 |
amult=weedthresh |
| 27 |
asmult=weedthresh |
| 28 |
|
| 29 |
if (icasea le 2) then begin |
| 30 |
iasym=0 |
| 31 |
nm=24 |
| 32 |
paroff=7 |
| 33 |
sigoff=5 |
| 34 |
aoff=12 |
| 35 |
endif else begin |
| 36 |
iasym=1 |
| 37 |
nm=26 |
| 38 |
paroff=8 |
| 39 |
sigoff=6 |
| 40 |
aoff=14 |
| 41 |
end |
| 42 |
|
| 43 |
if (keyword_set(prev)) then begin |
| 44 |
m1=readx(prev,nlin=nm1) |
| 45 |
m2=readx(infile,nlin=nm2) |
| 46 |
z=dblarr(nm,310,50) |
| 47 |
for i=0,nm1-1 do z(*,m1(0,i),m1(1,i))=m1(*,i) |
| 48 |
for i=0,nm2-1 do z(*,m2(0,i),m2(1,i))=m2(*,i) |
| 49 |
m3=reform(z,nm,15500) |
| 50 |
w=where(m3[2,*] ne 0,nm3) |
| 51 |
m3=m3[*,w] |
| 52 |
|
| 53 |
p3=getmodel(m3, icasea) |
| 54 |
lout=reform(m3[0,*]) |
| 55 |
nout=reform(m3[1,*]) |
| 56 |
ix=lonarr(max(lout)+1,max(nout)+1) |
| 57 |
for i=0,nm3-1 do ix(lout(i),nout(i))=i |
| 58 |
ix1=lonarr(nm2) |
| 59 |
for i=0,nm2-1 do ix1(i)=ix(m2(0,i),m2(1,i)) |
| 60 |
p2=p3[*,ix1] |
| 61 |
|
| 62 |
m=m2 |
| 63 |
p=p2 |
| 64 |
endif else begin |
| 65 |
m=readx(infile,nlin=nm2) |
| 66 |
p=getmodel(m, icasea) |
| 67 |
endelse |
| 68 |
|
| 69 |
gind=lonarr(nm2) |
| 70 |
csig=sqrt(binmics/m[4,*] > 1) |
| 71 |
wp=where(abs(p[2,*]-m[2,*]) lt numult*m[2+sigoff,*]*csig and $ |
| 72 |
abs(p[3,*]-m[3,*]) lt ampmult*m[3+sigoff,*]*csig and $ |
| 73 |
abs(p[4,*]-m[4,*]) lt widmult*m[4+sigoff,*]*csig and $ |
| 74 |
abs(1000*p[5,*]-m[aoff,*]) lt amult*m[aoff+6,*]*csig and $ |
| 75 |
abs(1000*p[6,*]-m[1+aoff,*]) lt amult*m[aoff+7,*]*csig and $ |
| 76 |
abs(1000*p[7,*]-m[2+aoff,*]) lt amult*m[aoff+8,*]*csig and $ |
| 77 |
abs(1000*p[8,*]-m[3+aoff,*]) lt amult*m[aoff+9,*]*csig and $ |
| 78 |
abs(1000*p[9,*]-m[4+aoff,*]) lt amult*m[aoff+10,*]*csig and $ |
| 79 |
abs(1000*p[10,*]-m[5+aoff,*]) lt amult*m[aoff+11,*]*csig) |
| 80 |
|
| 81 |
gind[wp]=1 |
| 82 |
if (iasym ne 0) then begin |
| 83 |
wa=where(abs(p[12,*]-m[7,*]) gt asmult*m[7+sigoff,*]*csig or $ |
| 84 |
m[7,*] lt 0,nwa) |
| 85 |
if (nwa gt 0) then gind[wa]=0 |
| 86 |
endif |
| 87 |
wg=where(gind,nwg,ncomp=nbad) |
| 88 |
|
| 89 |
print,'# modes rejected in iteration 1:',nbad |
| 90 |
if (iasym ne 0) then print,'# from asym:',nwa |
| 91 |
|
| 92 |
if (keyword_set(prev)) then begin |
| 93 |
m2=m[*,wg] |
| 94 |
nm2=nwg |
| 95 |
z=dblarr(nm,310,50) |
| 96 |
for i=0,nm1-1 do z(*,m1(0,i),m1(1,i))=m1(*,i) |
| 97 |
for i=0,nm2-1 do z(*,m2(0,i),m2(1,i))=m2(*,i) |
| 98 |
m3=reform(z,nm,15500) |
| 99 |
w=where(m3[2,*] ne 0,nm3) |
| 100 |
m3=m3[*,w] |
| 101 |
|
| 102 |
p3=getmodel(m3, icasea) |
| 103 |
lout=reform(m3[0,*]) |
| 104 |
nout=reform(m3[1,*]) |
| 105 |
ix=lonarr(max(lout)+1,max(nout)+1) |
| 106 |
for i=0,nm3-1 do ix(lout(i),nout(i))=i |
| 107 |
ix1=lonarr(nm2) |
| 108 |
for i=0,nm2-1 do ix1(i)=ix(m2(0,i),m2(1,i)) |
| 109 |
p2=p3[*,ix1] |
| 110 |
|
| 111 |
m=m2 |
| 112 |
p=p2 |
| 113 |
endif else begin |
| 114 |
m=m[*,wg] |
| 115 |
p=getmodel(m, icasea) |
| 116 |
endelse |
| 117 |
|
| 118 |
gind=lonarr(nwg) |
| 119 |
csig=sqrt(binmics/m[4,*] > 1) |
| 120 |
|
| 121 |
wp=where(abs(p[2,*]-m[2,*]) lt numult*m[2+sigoff,*]*csig and $ |
| 122 |
abs(p[3,*]-m[3,*]) lt ampmult*m[3+sigoff,*]*csig and $ |
| 123 |
abs(p[4,*]-m[4,*]) lt widmult*m[4+sigoff,*]*csig and $ |
| 124 |
abs(1000*p[5,*]-m[aoff,*]) lt amult*m[aoff+6,*]*csig and $ |
| 125 |
abs(1000*p[6,*]-m[1+aoff,*]) lt amult*m[aoff+7,*]*csig and $ |
| 126 |
abs(1000*p[7,*]-m[2+aoff,*]) lt amult*m[aoff+8,*]*csig and $ |
| 127 |
abs(1000*p[8,*]-m[3+aoff,*]) lt amult*m[aoff+9,*]*csig and $ |
| 128 |
abs(1000*p[9,*]-m[4+aoff,*]) lt amult*m[aoff+10,*]*csig and $ |
| 129 |
abs(1000*p[10,*]-m[5+aoff,*]) lt amult*m[aoff+11,*]*csig) |
| 130 |
|
| 131 |
gind[wp]=1 |
| 132 |
if (iasym ne 0) then begin |
| 133 |
wa=where(abs(p[12,*]-m[7,*]) gt asmult*m[7+sigoff,*]*csig or $ |
| 134 |
m[7,*] lt 0,nwa) |
| 135 |
if (nwa gt 0) then gind[wa]=0 |
| 136 |
endif |
| 137 |
wg=where(gind,nwg,ncomp=nbad) |
| 138 |
|
| 139 |
print,'# modes rejected in iteration 2:',nbad |
| 140 |
if (iasym ne 0) then print,'# from asym:',nwa |
| 141 |
|
| 142 |
m=m[*,wg] |
| 143 |
openw,1,outfile |
| 144 |
for i=0,nwg-1 do begin |
| 145 |
printf,1,m[*,i],format='(i3,i4,f11.4,99g16.8)' |
| 146 |
endfor |
| 147 |
close,1 |
| 148 |
|
| 149 |
END |
