nt=1 y6=dblarr(24,310,35,nt) y18=dblarr(48,310,35,nt) y36=dblarr(84,310,35,nt) par=dblarr(12,310,35,nt) a0=dblarr(9,310,35) names1=strarr(nt) names2=strarr(nt) nms6=lonarr(nt) nms18=lonarr(nt) nms36=lonarr(nt) openr,1,'list' for it=0,nt-1 do begin nm6=0l nm18=0l nm36=0l q='xxx' readf,1,nm6,nm18,nm36,q q1=str_sep(q,' ') name1=q1(1) name2=q1(2) print,it,nm6,nm18,nm36,' ',name1,' ',name2 names1(it)=name1 names2(it)=name2 nms6(it)=nm6 nms18(it)=nm18 nms36(it)=nm36 q=dblarr(24,nm6) openr,2,'m'+name1+'q.'+name2 readf,2,q close,2 for im=0,nm6-1 do y6(*,q(0,im),q(1,im),it)=q(*,im) q=dblarr(48,nm18) openr,2,'m'+name1+'q.'+name2+'.18' readf,2,q close,2 for im=0,nm18-1 do y18(*,q(0,im),q(1,im),it)=q(*,im) q=dblarr(84,nm36) openr,2,'m'+name1+'q.'+name2+'.36' readf,2,q close,2 for im=0,nm36-1 do y36(*,q(0,im),q(1,im),it)=q(*,im) q=dblarr(12,4724) openr,2,'par'+name1+'.'+name2 readf,2,q close,2 wp=where((q(1,*) le 300) and (q(0,*) le 34),np) for ip=0,np-1 do par(*,q(1,wp(ip)),q(0,wp(ip)),it)=q(*,wp(ip)) end close,1 ; Get splittings pridicted from 360d inversion q=dblarr(9,4627) openr,1,'a0' readf,1,q close,1 wa=where((q(0,*) le 300) and (q(1,*) le 34),na) for i=0,na-1 do a0(*,q(0,wa(i)),q(1,wa(i)))=q(*,wa(i)) nx=310l*35l*nt y6=reform(y6,24,nx,/overwrite) y18=reform(y18,48,nx,/overwrite) y36=reform(y36,84,nx,/overwrite) par=reform(par,12,nx,/overwrite) a0=reform(a0,9,310l*35l,/overwrite) ; Old code ;a0x=reform(rebin(reform(a0,9,310l*35l,1),9,310l*35l,nt),9,nx) a0x=reform(rebin(reform(a0,9,310l*35l,1),9,310l*35l,nt),9,310l,35l,nt) df=1e6/72./86400. l=reform(par(1,*)) n=reform(par(0,*)) f=reform(par(2,*)) mask0=intarr(nx) mask0(where(y6(2,*) ne 0))=1 na=6 mask6=intarr(nx) ; Toss modes more than 0.25 sigma away from input guess w6=where(abs((y6(2,*)-par(2,*))/y6(7,*)) le 0.25) mask6(w6)=1 ; Toss modes with large errors given width. Most likely low S/N. ; Added nwx test 020625 wx=where(y6(7,*)/y6(4,*)^0.5*(l+0.5)^0.5 ge 0.5,nwx) if (nwx gt 0) then mask6(wx)=0 ; Factor to jack up errors for modes with measured widths less than 1/72d. csig=reform(sqrt(0.16/y6(4,*)>1)) y6a=y6 y6a(7,*)=y6(7,*)*csig for ia=0,na-1 do y6a(12+na+ia,*)=y6(12+na+ia,*)*csig ; Toss modes more than 10 sigma away from prediction ;rx=dblarr(nx) ;for ia=0,na-2,2 do rx=rx>abs((y6a(12+ia,*)-a0x(3+ia,*))/y6a(12+na+ia,*)) ;mask6(where(rx gt 10))=0 ; New code y6a=reform(y6a,24,310l,35,nt) rx=dblarr(310l,35l,nt) for it=0,nt-1 do begin for ia=0,na-2,2 do begin w0=where(y6a(2,*,0,it) ne 0) mmm=median(y6a(12+ia,w0,0,it)-a0x(3+ia,w0,0,it)) q=reform(abs((y6a(12+ia,*,*,it)-a0x(3+ia,*,*,it)-mmm)/y6a(12+na+ia,*,*,it))) q(where(y6a(12+na+ia,*,*,it) eq 0))=0 rx(*,*,it)=rx(*,*,it)>q end end wq=where(rx gt 10,nq) if (nq gt 0) then mask6(wq)=0 ; Toss low freq isolated mode(s) mask6(where(f le 900))=0 w6=where(mask6 eq 1) na=18 mask18=mask6 mask18(where(y18(2,*) eq 0))=0 y18a=y18 y18a(7,*)=y18(7,*)*csig for ia=0,na-1 do y18a(12+na+ia,*)=y18(12+na+ia,*)*csig ; Toss modes with much larger errors or large changes from 6 a-coeff case. ry=reform(abs((y18(2,*)-y6(2,*))/y18a(7,*))) for ia=0,5 do ry=ry>reform(abs((y18(12+ia,*)-y6(12+ia,*))/y18a(12+na+ia))) sy=reform(y18(7,*)/y6(7,*)) for ia=0,5 do sy=sy>reform(y18(12+na+ia,*)/y6(18+ia,*)) mask18(where((ry ge 2) or (sy ge 2)))=0 w18=where(mask18 eq 1) na=36 mask36=mask18 mask36(where(y36(2,*) eq 0))=0 y36a=y36 y36a(7,*)=y36(7,*)*csig for ia=0,na-1 do y36a(12+na+ia,*)=y36(12+na+ia,*)*csig ; Toss modes with much larger errors or large changes from 6 a-coeff case. rz=reform(abs((y36(2,*)-y6(2,*))/y36a(7,*))) for ia=0,5 do rz=rz>reform(abs((y36(12+ia,*)-y6(12+ia,*))/y36a(12+na+ia))) sz=reform(y36(7,*)/y6(7,*)) for ia=0,5 do sz=sz>reform(y36(12+na+ia,*)/y6(18+ia,*)) mask36(where((rz ge 2) or (sz ge 2)))=0 w36=where(mask36 eq 1) ; Write things out na=6 y6a=reform(y6a,12+2*na,310l,35l,nt,/overwrite) mask6=reform(mask6,310l,35l,nt,/overwrite) for it=0,nt-1 do begin openw,1,'m'+names1(it)+'qr.'+names2(it) for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask6(il,in,it) eq 1) then printf,1,y6a(*,il,in,it),format='(2i4,f10.4,f10.4,4f8.4,2f9.4,f10.4,f9.4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'freq'+names1(it)+'qr.'+names2(it) for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask6(il,in,it) eq 1) then printf,1,y6a([0,1,2,7],il,in,it),format='(2i4,f10.4,f8.4)' end end close,1 end i1=[0,1,2,12+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it) for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask6(il,in,it) eq 1) then printf,1,y6a(i1,il,in,it),format='(2i4,f10.4,f10.4,99g14.6)' end end close,1 end i1=[0,1,2,13+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'splite'+names1(it)+'qr.'+names2(it) for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask6(il,in,it) eq 1) then printf,1,y6a(i1,il,in,it),format='(2i4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.2d' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask6(il,in,it) eq 1) then begin for ia=0,na-2,2 do begin if y6a(12+na+ia,il,in,it) ne 0 then begin printf,1,y6a([0,1,2],il,in,it),ia/2+1,3,na/2,y6a([12+ia,12+na+ia],il,in,it),format='(i4,i3,f8.2,i4,i2,i3,2g14.6)' end end end end end close,1 end na=18 y18a=reform(y18a,12+2*na,310l,35l,nt,/overwrite) mask18=reform(mask18,310l,35l,nt,/overwrite) for it=0,nt-1 do begin openw,1,'m'+names1(it)+'qr.'+names2(it)+'.18' for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask18(il,in,it) eq 1) then printf,1,y18a(*,il,in,it),format='(2i4,f10.4,f10.4,4f8.4,2f9.4,f10.4,f9.4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'freq'+names1(it)+'qr.'+names2(it)+'.18' for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask18(il,in,it) eq 1) then printf,1,y18a([0,1,2,7],il,in,it),format='(2i4,f10.4,f8.4)' end end close,1 end i1=[0,1,2,12+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.18' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask18(il,in,it) eq 1) then printf,1,y18a(i1,il,in,it),format='(2i4,f10.4,f10.4,99g14.6)' end end close,1 end i1=[0,1,2,13+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'splite'+names1(it)+'qr.'+names2(it)+'.18' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask18(il,in,it) eq 1) then printf,1,y18a(i1,il,in,it),format='(2i4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.18.2d' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask18(il,in,it) eq 1) then begin for ia=0,na-2,2 do begin if y18a(12+na+ia,il,in,it) ne 0 then begin printf,1,y18a([0,1,2],il,in,it),ia/2+1,3,na/2,y18a([12+ia,12+na+ia],il,in,it),format='(i4,i3,f8.2,i4,i2,i3,2g14.6)' end end end end end close,1 end na=36 y36a=reform(y36a,12+2*na,310l,35l,nt,/overwrite) mask36=reform(mask36,310l,35l,nt,/overwrite) for it=0,nt-1 do begin openw,1,'m'+names1(it)+'qr.'+names2(it)+'.36' for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask36(il,in,it) eq 1) then printf,1,y36a(*,il,in,it),format='(2i4,f10.4,f10.4,4f8.4,2f9.4,f10.4,f9.4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'freq'+names1(it)+'qr.'+names2(it)+'.36' for il=0l,310l-1 do begin for in=0l,35l-1 do begin if (mask36(il,in,it) eq 1) then printf,1,y36a([0,1,2,7],il,in,it),format='(2i4,f10.4,f8.4)' end end close,1 end i1=[0,1,2,12+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.36' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask36(il,in,it) eq 1) then printf,1,y36a(i1,il,in,it),format='(2i4,f10.4,f10.4,99g14.6)' end end close,1 end i1=[0,1,2,12+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.36f' for il=1l,310l-1 do begin for in=0l,0l do begin if (mask36(il,in,it) eq 1) then printf,1,y36a(i1,il,in,it),format='(2i4,f10.4,f10.4,99g14.6)' end end close,1 end i1=[0,1,2,13+2*lindgen(na)] for it=0,nt-1 do begin openw,1,'splite'+names1(it)+'qr.'+names2(it)+'.36' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask36(il,in,it) eq 1) then printf,1,y36a(i1,il,in,it),format='(2i4,f10.4,99g14.6)' end end close,1 end for it=0,nt-1 do begin openw,1,'split'+names1(it)+'qr.'+names2(it)+'.36.2d' for il=1l,310l-1 do begin for in=0l,35l-1 do begin if (mask36(il,in,it) eq 1) then begin for ia=0,na-2,2 do begin if y36a(12+na+ia,il,in,it) ne 0 then begin printf,1,y36a([0,1,2],il,in,it),ia/2+1,3,na/2,y36a([12+ia,12+na+ia],il,in,it),format='(i4,i3,f8.2,i4,i2,i3,2g14.6)' end end end end end close,1 end end