| 1 |
PRO MKSEC, firstday, ndt, tstep, ref=ref |
| 2 |
|
| 3 |
; this code assumes all timeseries begin and end on day boundaries. |
| 4 |
; also assumes all timeseries use same t_start_epoch and that |
| 5 |
; t_start_step is specified in seconds, but t_start_step can vary. |
| 6 |
|
| 7 |
if (n_elements(ref) eq 0) then ref='mdi.secs_72d' |
| 8 |
spawn,'show_info -q '+ref+ '\[] key=t_start_index,t_start_step,ndt,t_step > daylist.tmp' |
| 9 |
|
| 10 |
ndays=long(ndt*tstep/86400) ;number of days in requested timeseries |
| 11 |
lastday=firstday+ndays-1 ;last day of requested timeseries |
| 12 |
ppd=long(86400/tstep) ;points per day |
| 13 |
|
| 14 |
|
| 15 |
daystr=strtrim(string(firstday),2) |
| 16 |
ndstr=strtrim(string(ndays),2) |
| 17 |
file='secs.'+ndstr+'.'+daystr |
| 18 |
|
| 19 |
arr=readx('daylist.tmp',nl=nl) |
| 20 |
test=reform(arr[3,*]) |
| 21 |
wtest=where(test ne tstep,nwt) |
| 22 |
if (nwt gt 0) then begin |
| 23 |
print,'ERROR: reference timeseries has wrong t_step' |
| 24 |
return |
| 25 |
endif |
| 26 |
file_delete,'daylist.tmp' |
| 27 |
|
| 28 |
fd=long(arr[0,*]) |
| 29 |
t_start_step=long(arr[1,0]) |
| 30 |
fd=fd[*]*t_start_step/86400 |
| 31 |
nd=long(arr[2,*]) |
| 32 |
nd=nd[*]*tstep/86400 |
| 33 |
ld=fd+nd-1 |
| 34 |
;fd contains the first day of each refence timeseries, nd the number of days in each, |
| 35 |
;ld the last day of each. |
| 36 |
|
| 37 |
;first check if requested section file is identical to one of the reference section files |
| 38 |
wt=where(fd eq firstday and nd eq ndays,nwt) |
| 39 |
if (nwt eq 1) then begin |
| 40 |
spawn,'show_info -q '+ref+ '\['+strtrim(string(firstday),2)+'d] key=nsecs > '+file |
| 41 |
spawn,'show_info -q '+ref+ '\['+strtrim(string(firstday),2)+'d] key=secs >> '+file |
| 42 |
return |
| 43 |
endif |
| 44 |
|
| 45 |
;check for no overlap between requested section file and reference section files |
| 46 |
if (firstday gt ld[nl-1] or lastday lt fd[0]) then begin |
| 47 |
print,"ERROR: no overlap between requested section file and reference section files" |
| 48 |
return |
| 49 |
endif |
| 50 |
|
| 51 |
w1=where(firstday ge fd and firstday le ld,nw1) |
| 52 |
if (nw1 eq 0) then w1=(where(firstday lt fd))[0] |
| 53 |
w2=where(lastday ge fd and lastday le ld,nw2) |
| 54 |
if (nw2 eq 0) then w2=(where(lastday gt ld,nw))[nw-1] |
| 55 |
|
| 56 |
;in this case the requested section file lies entirely between 2 reference section files |
| 57 |
if (w2 lt w1) then return |
| 58 |
|
| 59 |
ts=intarr(ndt) |
| 60 |
nsecs=0L |
| 61 |
rec=lonarr(2) |
| 62 |
;w1 and w2 contain 1 element each, use total() to convert them from arrays to numbers |
| 63 |
for i=total(w1),total(w2) do begin |
| 64 |
sub=intarr(nd[i]*ppd) |
| 65 |
spawn,'show_info -q '+ref+ '\['+strtrim(string(fd[i]),2)+'d] key=nsecs',hold |
| 66 |
reads,hold,nsecs |
| 67 |
spawn,'show_info -q '+ref+ '\['+strtrim(string(fd[i]),2)+'d] key=secs',hold |
| 68 |
for j=0,nsecs-1 do begin |
| 69 |
reads,hold[j],rec |
| 70 |
sub[rec[0]:rec[1]]=1 |
| 71 |
endfor |
| 72 |
|
| 73 |
if (firstday le fd[i]) then begin |
| 74 |
tsind1=(fd[i] - firstday)*ppd |
| 75 |
sbind1=0 |
| 76 |
endif else begin |
| 77 |
tsind1=0 |
| 78 |
sbind1=(firstday - fd[i])*ppd |
| 79 |
endelse |
| 80 |
if (lastday ge ld[i]) then begin |
| 81 |
tsind2=ndt - (lastday - ld[i])*ppd - 1 |
| 82 |
sbind2=nd[i]*ppd - 1 |
| 83 |
endif else begin |
| 84 |
tsind2=ndt - 1 |
| 85 |
sbind2=nd[i]*ppd - (ld[i] - lastday)*ppd - 1 |
| 86 |
endelse |
| 87 |
|
| 88 |
ts[tsind1:tsind2]=sub[sbind1:sbind2] |
| 89 |
endfor |
| 90 |
|
| 91 |
ws=where(ts) |
| 92 |
as=arrsumm(ws) |
| 93 |
; the following line removes sections of less than 100 time points |
| 94 |
index=where((as[1,*]-as[0,*]) gt 100,nsecs) |
| 95 |
secs=as[*,index] |
| 96 |
|
| 97 |
openw,1,file |
| 98 |
printf,1,nsecs |
| 99 |
printf,1,secs |
| 100 |
close,1 |
| 101 |
|
| 102 |
END |
