c . 1 . 2 . 3 . 4 . 5 . 6 . 7.. c c linear combination of six signals, to fit a seventh signal c May be used for the determination of a moment tensor c Oct. 2009 EW c program linreg6 use nas_system parameter(ndim=16500) dimension r(ndim),s(ndim),t(ndim),syn(ndim) dimension x(ndim),y(ndim),z(ndim),w(ndim), & wr(ndim),aik(6,6),rs(6),f(6) character*48 inf1,inf2,inf3,inf4,inf5,inf6,inf7,synf,resf,ddir character*20 form character*72 text1,text2,text3,text4,text5,text6,text7 write(6,2) write(6,2) ' LINREG6: solve w = a*r +b*s +c*t +d*x +e*y +f*z' write(6,2) ' for a, b, c, d, e,f ' write(6,2) ' r=file1, s = file2, etc.' write(6,2) ' the output file is the residue w -a*r -b*s -...' write(6,2) ' The parameterfile linreg.par must contain' write(6,2) ' name of the data directory' write(6,2) ' names of the r, s, t, x, y, z files, one line each' write(6,2) ' name of the w file' write(6,2) ' eps: stabilization weight, can normally be zero' write(6,2) ' mn: remove mean if =1, remove trend if =2' 2 format(a) open(11,file='linreg.par',status='old') read(11,*) ddir read(11,*) inf1 read(11,*) inf2 read(11,*) inf3 read(11,*) inf4 read(11,*) inf5 read(11,*) inf6 read(11,*) inf7 read(11,*) eps read(11,*) mn n=ndim inf1=trim(ddir)//'\'//trim(inf1) inf2=trim(ddir)//'\'//trim(inf2) inf3=trim(ddir)//'\'//trim(inf3) inf4=trim(ddir)//'\'//trim(inf4) inf5=trim(ddir)//'\'//trim(inf5) inf6=trim(ddir)//'\'//trim(inf6) inf7=trim(ddir)//'\'//trim(inf7) c read data call input(inf1,text1,n,form,dt,tmin,tsec,r) call input(inf2,text2,n,form,dt,tmin,tsec,s) call input(inf3,text4,n,form,dt,tmin,tsec,t) call input(inf4,text4,n,form,dt,tmin,tsec,x) call input(inf5,text5,n,form,dt,tmin,tsec,y) call input(inf6,text6,n,form,dt,tmin,tsec,z) call input(inf7,text7,n,form,dt,tmin,tsec,w) if(mn.eq.1) then call mean(r,n) call mean(s,n) call mean(t,n) call mean(x,n) call mean(y,n) call mean(z,n) call mean(w,n) else if (mn.eq.2) then call trend(r,n) call trend(s,n) call trend(t,n) call trend(x,n) call trend(y,n) call trend(z,n) call trend(w,n) endif do j=1,6 rs(j)=0. do k=j,6 aik(j,k)=0. enddo aik(j,j)=eps enddo do i=1,n aik(1,1)=aik(1,1)+r(i)*r(i) aik(1,2)=aik(1,2)+r(i)*s(i) aik(1,3)=aik(1,3)+r(i)*t(i) aik(1,4)=aik(1,4)+r(i)*x(i) aik(1,5)=aik(1,5)+r(i)*y(i) aik(1,6)=aik(1,6)+r(i)*z(i) aik(2,2)=aik(2,2)+s(i)*s(i) aik(2,3)=aik(2,3)+s(i)*t(i) aik(2,4)=aik(2,4)+s(i)*x(i) aik(2,5)=aik(2,5)+s(i)*y(i) aik(2,6)=aik(2,6)+s(i)*z(i) aik(3,3)=aik(3,3)+t(i)*t(i) aik(3,4)=aik(3,4)+t(i)*x(i) aik(3,5)=aik(3,5)+t(i)*y(i) aik(3,6)=aik(3,6)+t(i)*z(i) aik(4,4)=aik(4,4)+x(i)*x(i) aik(4,5)=aik(4,5)+x(i)*y(i) aik(4,6)=aik(4,6)+x(i)*z(i) aik(5,5)=aik(5,5)+y(i)*y(i) aik(5,6)=aik(5,6)+y(i)*z(i) aik(6,6)=aik(6,6)+z(i)*z(i) rs(1)= rs(1)+r(i)*w(i) rs(2)= rs(2)+s(i)*w(i) rs(3)= rs(3)+t(i)*w(i) rs(4)= rs(4)+x(i)*w(i) rs(5)= rs(5)+y(i)*w(i) rs(6)= rs(6)+z(i)*w(i) enddo do i=1,5 do k=i+1,6 aik(k,i)=aik(i,k) enddo enddo call gauss(aik,6,6,rs,f) wq=0. wrq=0. do i=1,n syn(i)=f(1)*r(i)+f(2)*s(i)+f(3)*t(i)+f(4)*x(i)+f(5)*y(i)+f(6)*z(i) wr(i)=w(i)-syn(i) wq = wq+ w(i)**2 wrq=wrq+wr(i)**2 enddo write(6,1) trim(inf7),trim(inf1),f(1),trim(inf2),f(2), & trim(inf3),f(3),trim(inf4),f(4),trim(inf5),f(5),trim(inf6), & f(6),wq,wrq,100.*sqrt(wrq/wq) 1 format(/' Decomposition of signal w = ',a,':'/ & ' coefficient of series r = ',a,f13.6/ & ' coefficient of series s = ',a,f13.6/ & ' coefficient of series t = ',a,f13.6/ & ' coefficient of series x = ',a,f13.6/ & ' coefficient of series y = ',a,f13.6/ & ' coefficient of series z = ',a,f13.6// & ' Energy of series w ',e10.3/ & ' Energy of residual ',e10.3/ & ' rms residual ',f10.6,' %') synf=trim(inf7)//'.syn' call output(synf,text1,n,form,dt,tmin,tsec,syn) resf=trim(inf7)//'.res' call output(resf,text1,n,form,dt,tmin,tsec,wr) c compose winplot.par file open(12,file='winplot.par') write(12,*) '0, 9, 24., 16., 1, 0, 0.8' write(12,'(a)') trim(inf7) write(12,'(a)') trim(synf) write(12,'(a)') trim(resf) write(12,'(a)') trim(inf1) write(12,'(a)') trim(inf2) write(12,'(a)') trim(inf3) write(12,'(a)') trim(inf4) write(12,'(a)') trim(inf5) write(12,'(a)') trim(inf6) write(12,'(a)') 'linreg6 synthesis' close(12) stop end c subroutine mean(x,n) c remove average dimension x(n) sum=0. do 1 j=1,n 1 sum=sum+x(j) sum=sum/n do 2 j=1,n 2 x(j)=x(j)-sum return end c subroutine trend(x,n) c remove linear trend dimension x(n) gn = float(n) alpha = 0.5*gn*(gn+1.) beta = (2.*gn+1.)*(gn+1.)*gn/6. det = gn*beta-alpha*alpha sx = 0. sjx = 0. do 1 j=1,n sx = sx+x(j) 1 sjx = sjx+x(j)*float(j) a = (sx*beta-sjx*alpha)/det b = (sjx*gn-sx*alpha)/det do 2 j=1,n 2 x(j) = x(j)-a-b*float(j) return end c subroutine input(name,text,n,form,dt,tmin,tsec,x) dimension x(n) character name*(*),form*20,text*72 call inpasc(name,text,n,form,dt,tmin,tsec,x) return end c subroutine inpasc(name,text,n,form,dt,tmin,tsec,x) dimension x(n) character form*20,name*(*),text*72,pro*60 open(7,file=name,status='old') read(7,'(a)') text 21 read(7,'(a)') pro if(pro(1:1).eq.'%') goto 21 read(pro,1) nn,form,dt,tmin,tsec 1 format(i10,a20,3f10.3) write(6,22) trim(name),nn,n 22 format(' file ',a,' wird gelesen:',i7,' /',i7,' punkte') if(n.lt.nn.or.nn.eq.0) then c write(6,'("punktzahl auf",i7," festgesetzt")') n nn=n endif n=nn if(form(1:4).eq.'frei') then read(7,*,err=25,end=23) (x(j),j=1,n) else read(7,form,end=23) (x(j),j=1,n) endif close(7) return 23 n=j-1 write(6,24) n 24 format('Punktzahl ',i8,' durch Datenende bestimmt') close(7) return 25 write(6,26) name,j 26 format(' Input error reading file ',a,' at sample # ',i7) stop end c subroutine output(name,text,n,form,dt,tmin,tsec,x) dimension x(n) character name*(*),text*72,form*20 xmax=0. do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) nvor=int(log10(max(xmax,2.))+1.) ndec=max(0,10-nvor) i1=index(form,'(') i2=index(form,')') if(index(form,'frei').gt.0) then write(form,'(a,a,a)') '(5f13.',ndec,')' else if(i1.lt.1.or.i2.lt.i1+5) then write(6,'(a,a,a)') ' Ausgabeformat ',form, &' nicht interpretierbar' stop endif write(6,5) trim(name) 5 format(/' writing residue into file ',a) open(8,file=name) write(8,'(a)') text write(8,1) n,form,dt,tmin,tsec 1 format(i10,a20,f10.3,2f10.3) write(8,form) (x(j),j=1,n) close(8) return end subroutine gauss(aik,m,n,rs,f) c solve linear equations dimension aik(n,n),rs(n),f(n),h(14),imax(13) do 1401 j=1,m aikmax=0. do 1402 k=1,m h(k)=aik(j,k) if(abs(h(k)).le.aikmax) go to 1402 aikmax=abs(h(k)) index=k 1402 continue h(m+1)=rs(j) do 1403 k=1,m q=aik(k,index)/h(index) do 1404 l=1,m 1404 aik(k,l)=aik(k,l)-q*h(l) 1403 rs(k)=rs(k)-q*h(m+1) do 1405 k=1,m 1405 aik(j,k)=h(k) rs(j)=h(m+1) 1401 imax(j)=index do 1406 j=1,m index=imax(j) 1406 f(index)=rs(j)/aik(j,index) return end