c . 1 . 2 . 3 . 4 . 5 . 6 . 7.. c c linear combination of two signals, to fit a third one c July 1996 e.w. c program lincomb use nas_system parameter(ndim=60000) dimension x(ndim),y(ndim),z(ndim),zz(ndim) character inf1*24,inf2*24,inf3*24,outf*28,form*20 character*72 text1,text2,text3 write(6,*) write(6,'(a,a)') 'lincomb2 file1 file2 file3 [n1 n2]:', & ' solve z = a*x + b*y + residue' write(6,'(a,a)') 'where x=signal1, y=signal2, z=signal3, and the', & ' constants a and b are ' write(6,'(a,a)') 'determined for best fit. A window may be speci', & 'fied with n1 and n2.' write(6,'(a)') 'The residue is written into file3.res' narg=iargc() if(narg.lt.3) stop call getarg(1,inf1) call getarg(2,inf2) call getarg(3,inf3) if (narg.eq.5) then call getarg(4,form) read(form,*) nstart call getarg(5,form) read(form,*) nstop else nstart=1 nstop=60000 endif c read data n1=ndim call input(inf1,text1,n1,form,dt,tmin,tsec,x) n2=ndim call input(inf2,text2,n2,form,dt,tmin,tsec,y) n3=ndim call input(inf3,text3,n3,form,dt,tmin,tsec,z) n=min(n1,n2,n3) nstart=max(nstart,1) nstop=min(nstop,n) call mean(x,n) call mean(y,n) call mean(z,n) a1=0. a2=0. a3=0. b1=0. b2=0. do i=nstart,nstop a1=a1+x(i)**2 a2=a2+y(i)**2 a3=a3+x(i)*y(i) b1=b1+x(i)*z(i) b2=b2+y(i)*z(i) enddo det=a1*a2-a3*a3 a=(b1*a2-b2*a3)/det b=(b2*a1-b1*a3)/det do i=1,n zz(i)=z(i)-a*x(i)-b*y(i) enddo call mean(zz,n) xq=0. yq=0. zq=0. zzq=0. do i=nstart,nstop xq=xq+x(i)**2 yq=yq+y(i)**2 zq=zq+z(i)**2 zzq=zzq+zz(i)**2 enddo nf=nstop-nstart+1 xq=sqrt(xq/nf) yq=sqrt(yq/nf) zq=sqrt(zq/nf) zzq=sqrt(zzq/nf) write(6,1) a,b,xq,yq,zq,zzq,100.*zzq/zq,nstart,nstop 1 format(' Koeff. series x = ',f15.6,' series y = ',f15.6/ & ' rms Amp. series x = ',f15.6,' series y = ',f15.6/ & ' series z = ',f15.6,' residue = ',f15.6/ & ' the residue is ',f6.1,' % of the rms amplitude'/ & ' data were fitted from samples ',i7,' to ',i7, & ' (mean removed)') c Restsignal in File inf3.res schreiben text3='Rest' outf=trim(inf3)//'.res' call output(outf,text3,n,form,dt,tmin,tsec,zz) 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 form*20,name*24,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('reading file ',a,', ',i7,' /',i7,' samples') 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(' end of file encountered after ',i7,' samples') close(7) return 25 write(6,26) trim(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*28,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)') ' illegal format specification ',form stop endif write(6,5) trim(name),n,form 5 format('writing file ',a,i10,' samples in format ',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