PROGRAM TRICAL C This is a special version of CALEX involving a transformation of the C XYZ output signals to UVW before the analysis, and a transformation C of the results back to XYZ. See the program description and section C 5.6.7 of the new MSOP. C System calibration using arbitrary test signals C Simulation using impulse-invariant recursive filters C Reference: H.W.Schuessler, A signalprocessing approach to simulation, C Frequenz 35 (1981), 174-184 c This software is only applicable to calibration procedures where both c the stimulus and the response are digitally recorded. Since normally c both signals will be recorded with identical digitizers, the digitizer c response is irrelevant and will not be modelled. c According to Schuessler, the method can only be applied to systems c whose passband lies entirely within the Nyquist bandwidth. That is, c the impulse response of the system must have negligible energy outside c this bandwidth. Seismometers do in general not satisfy this condition, c especially when only the lower corner frequency is modelled so that c the seismometer appears as a high-pass filter. We therefore limit the c bandwidth of the system numerically with a recursive filter whose c corner period is specified by the alias parameter, and model the c real system plus the filter. Since the filter is known, its parameters c do not appear in the inversion as active parameters, and the user of c this program may not notice its presence except for a slight smoothing c of the synthetic output signal. To maintain visual consistency, the c plottable input and output signals, einf and ausf, are also filtered c with the anti-alias filter; however the inversion uses the original, c unfiltered input signal. c It is strongly recommended to plot and inspect the signals einf, ausf, c synt and rest after each run. synt is the synthetic output signal and c rest is the misfit. Under MS-Windows, the plot routine WINPLOT will c automatically produce a suitable plot. For details, see the program c description calex.doc . c The meaning of the 'amp' parameter is different for different experiments c (electrical or mechanical input, etc.). If a broadband-velocity seismometer c with the asymptotic (practically, mid-band) generator constant G is tested c on a shake table that has a displacement transducer, two possible cases are: c a) the displacement transducer has an electric output of D volts per meter c and its output signal is recorded with the same digitizer as the seismo- c meter output. Then amp = G / D. c b) The displacement transducer is digital, or has a separate digitizer, with c a responsivity A counts per meter. The seismometer output is digitized c with C counts per volt. Then amp = C * G / A. c Revisions: c 04 Dec 1998 c partial systems defined by keywords hp1, lp1, hp2, bp2, lp2 c anti-alias-filtered 'ausg' signal written into file 'ausf' c output in free format c 22 Feb 1999 c delay parameter alternatively renamed 'sub' and (ab)used for c compensation of direct input-to-output coupling (for calibration c of single-coil geophones in a half bridge). c c 06 Mar 2002 c the delay parameter is alternatively used for the compensation of tilt c effects on a shake table. c c 12 Mar 2002 c In place of using the delay parameter alternatively for the compensation c of galvanic coupling or tilt, the program now expects to find four c parameters with fixed meaning: amp for the gain, del for the delay, sub c for galvanic coupling, and til for the tilt. The sub parameter causes the c fraction sub of the input signal to be subtracted from the output signal. c The til parameter causes the twice-integrated signal to be added with c amplitude til (actually this is done with the output signal, but could be c done with the input signal as well, which would correspond better to the c physical situation on a shake table). c The parameter file is no longer compatible with that of the original c version. It must contain two additional lines with 'sub' and 'til' after c the 'del' line. c c til is in units of microradians per millimeter of table motion provided c that the motion of the table, the gravitational force coupled in by tilt, c and the sensitive axis of the seismometer have the same direction. More c generally, til is (tilt along the sensitive axis) per (table displacement c along the sensitive axis). c c 18 Sep 05 c names of input and output files are read from calex.par c The .par file must contain two additional lines with the file names c before all other parameters, and is no longer compatible with c earlier versions. File names should be enclosed in simple quotes c (apostrophs) especially when they contain spaces or the line contains c comments after the file names. Lines starting with a space are ignored. c c May 2006: ASL format is recognized. (Unreliably, as I found out later.) c July 2006: Name of parameter file can be entered as runstring. c Output file named after parameter file. c c August 2007: poles and zeros are listed after the inversion. c c October 2007: delay can be up to +-24*dt. The synthetic signal is computed c from sample ns1-24 to sample ns2+24 to allow for the time shift. DEL and c SUB can no longer be used at the same time. c c October 2008: the old "spaghetti" version of subroutine mini was replaced c by an equivalent but more comprehensible code. c c March - May 2009: several minor improvements. ns1 and ns2 are written into c the plot-parameter file so that only the modelled part of the signal is c plotted. The output files contain howewer the whole signal. c c March 2010: Reading ASL format did not work any more because the unreliable c FORTRAN free-format read statement was used. It was replaced by the reading c routine from SEIFE. Other programs had been updated earlier. c c ********** Version for homogeneous-triaxial sensors ********** March 2011 c reads aus3 output files, transforms the signals to UVW, applies CALEX to each c one, finally transforms parameters back to XYZ (Method of Peter Davis) c File names 'calex.*' chenged to 'trical.*' c c April 2015: start values of the normalized active parameters x were set to c zero only for the first component, are now zeroed for all three components. c (Thanks to Thomas Forbriger who noticed this normally inconsequential bug.) use nas_system c the 'use' statement is compiler specific, and not required for Linux parameter (mpar=12,msys=25,ns=480000) c needs heap 40000 kBytes for ns=480000, 60000 for ns=800000 implicit double precision (a-h,o-z) character endtxt*9, titel*72,name*3,typ*3,plop*12 character*50 eing,ausx,ausy,ausz character*5 einfile,ausf(3),synt(3),rest(3) character*1 uvw(3),xyz(3) dimension x(mpar) dimension ein(ns),aus(ns),aus3(ns,3),sum(ns),einf(ns),sta(ns) dimension csp(0:msys,3),cpp(msys,3,3) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m common /par/ x0(msys),rho(msys),x00(mpar),r00(mpar),typ(msys), & name(mpar) common /sys/ syspar(0:msys),sr(msys),si(msys),resid(msys), & residi(msys),m0,m1,m2 data zero,one/0.d0,1.d0/ data xyz,uvw /'X','Y','Z','U','V','W'/ c filenames of plottable signals einfile='einf3' ausf(1)='ausfu' ausf(2)='ausfv' ausf(3)='ausfw' synt(1)='syntu' synt(2)='syntv' synt(3)='syntw' rest(1)='restu' rest(2)='restv' rest(3)='restw' plop='winplot0.par' open(3,file='trical.par',status='old') c read header line and file names read(3,'(a)') titel open(4,file='trical.out') write(4,'(/,a,/)') titel write(6,'(/,a,/)') titel 201 read(3,*) eing ! only filenames are read here ! if(eing(1:1).eq.' '.or.eing(1:1).eq.'') goto 201 222 read(3,*) ausx if(ausx(1:1).eq.' '.or.ausx(1:1).eq.'') goto 222 223 read(3,*) ausy if(ausy(1:1).eq.' '.or.ausy(1:1).eq.'') goto 223 224 read(3,*) ausz if(ausz(1:1).eq.' '.or.ausz(1:1).eq.'') goto 224 close(3) C read data, transform and prefilter output signals c array 'ein' is used as temporary storage only n=ns call input(ausx,n,dt,ein,tma,tsa) do i=1,n aus3(i,1)=ein(i) enddo call input(ausy,n,dt,ein,tma,tsa) do i=1,n aus3(i,2)=ein(i) enddo call input(ausz,n,dt,ein,tma,tsa) do i=1,n aus3(i,3)=ein(i) enddo call trans(aus3,n) c ********************************************************************** c Loop over U, V, W output signals do 225 icomp=1,3 write(4,*) write(4,*) '***** analysing the ',uvw(icomp),' component *****' write(4,*) write(6,*) write(6,*) '***** analysing the ',uvw(icomp),' component *****' write(6,*) init=-1 mq=0 endtxt ='converged' do i=1,mpar x(i)=0.d0 enddo c corner period and order of the anti-alias filter call inpar(m,maxit,qac,finac,init,dt,ns1,ns2) ns1=1 ns2=n do i=1,n ein(i)=aus3(i,icomp) enddo call quadr(x,q,ein,aus,sum,einf,sta) c note: here (init<0) sum is the anti-alias filtered output signal and aus is ignored. do 28 i=1,n 28 aus(i)=sum(i) c now aus is the anti-alias filtered output signal call output(ausf(icomp),'ausg with anti-alias',n,dt,aus,tma,tsa) ne=ns call input(eing,ne,dtein,ein,tme,tse) ein0=ein(1) do 30 i=1,n 30 ein(i)=ein(i)-ein0 call quadr(x,qq,ein,aus,sum,einf,sta) do 31 i=1,n 31 einf(i)=sum(i) c einf is the anti-alias filtered input signal, used only for plotting if(icomp.eq.1) then call output(einfile,'eing with anti-alias',n,dt,einf,tme,tse) endif C prepare iteration, read start parameters iter=0 step=one noimp=0 init=0 call inpar(m,maxit,qac,finac,init,dt,ns1,ns2) if(ns1.le.0) ns1=1 if(ns2.le.0) ns2=n ns1=max(ns1,25) ns2=min(ns2,n-23) if(ns2.le.ns1) then write(6,'(a)') ' incorrect time window - check ns1, ns2' stop endif write(6,101) ns1,ns2 write(4,101) ns1,ns2 101 format(' data will be fitted from sample',i7,' to',i7) c start model qnorm=zero do 4 j=ns1,ns2 4 qnorm=qnorm+aus(j)**2 qn=qnorm call quadr(x,q,ein,aus,sum,einf,sta) c from here on (init=0) sum is the misfit do 96 i=1,n 96 sta(i)=zero do 97 i=ns1,ns2 97 sta(i)=aus(i)-sum(i) write(6,98) write(4,98) 98 format(' synthetics from start model:') call output(synt(icomp),'synth. for startmod.',n,dt,sta,tme,tse) write(6,99) write(4,99) 99 format(' misfit from start model:') call output(rest(icomp),'misfit for startmod.',n,dt,sum,tme,tse) write(6,103) (name(k),k=1,m) write(4,103) (name(k),k=1,m) 103 format(/' iter',9x,'RMS',5(9x,a3)/(17x,5(9x,a3))) write(6,104) write(4,104) if(mod(m,5).eq.0) then write(6,106) iter,dsqrt(q),(x00(k),k=1,m) write(4,106) iter,dsqrt(q),(x00(k),k=1,m) else write(6,104) iter,dsqrt(q),(x00(k),k=1,m) write(4,104) iter,dsqrt(q),(x00(k),k=1,m) endif 104 format(i5,6(1x,f11.6)/(17x,5(1x,f11.6))) 106 format(i5,6(1x,f11.6)) write(6,105) (r00(k),k=1,m) write(4,105) (r00(k),k=1,m) 105 format((15x,'+-',5(1x,f11.6))) write(6,104) write(4,104) C iteration with the conjugate-gradient method do 2 iter=1,maxit qalt=q call congrd(x,q,gnorm,iter,ein,aus,sum,einf,sta) if(mod(m,5).eq.0) then write(6,106) iter,dsqrt(q),(x(k),k=1,m) write(4,106) iter,dsqrt(q),(x(k),k=1,m) else write(6,104) iter,dsqrt(q),(x(k),k=1,m) write(4,104) iter,dsqrt(q),(x(k),k=1,m) endif c the iteration stops when in m steps, the fit was improved by less c than qac and the parameter vector was changed by less than finac. c The rms error is relative to the rms amplitude of the output signal. c Parameter corrections are relative to the search interval, rho. noimp=noimp+1 if((dsqrt(qalt)-dsqrt(q)).gt.qac.or.axi.gt.finac) noimp=0 if(noimp.ge.m) goto 5 2 continue endtxt='not conv.' iter=maxit 5 write(6,104) write(4,104) write(6,104) write(4,104) write(6,109) write(4,109) 109 format(' final system parameters:') do k=1,m cpp(k,icomp,1)=x00(k)+r00(k)*x(k) enddo write(6,103) (name(k),k=1,m) write(4,103) (name(k),k=1,m) write(6,104) iter,dsqrt(q),(x00(k)+r00(k)*x(k),k=1,m) write(4,104) iter,dsqrt(q),(x00(k)+r00(k)*x(k),k=1,m) write(6,108) mq,endtxt write(4,108) mq,endtxt 108 format(/' QUAD called',i5,' times: ',a) if(iter.lt.maxit) then do 196 i=1,n 196 sta(i)=zero do 197 i=ns1,ns2 197 sta(i)=aus(i)-sum(i) write(6,'(a)') ' output signal for final model...' write(4,'(a)') ' output signal for final model...' call output(synt(icomp),'synth. for final mod',n,dt,sta,tme,tse) endif do 6 i=1,n 6 aus(i)=zero do 7 i=ns1,ns2 7 aus(i)=sum(i) write(6,'(a)') ' residual error...' write(4,'(a)') ' residual error...' call output(rest(icomp),'obs.-synth.'//endtxt,n,dt,aus,tme,tse) c save poles for backtransformation csp(0,icomp)=syspar(0) do k=1,msys csp(k,icomp)=syspar(k) cpp(k,icomp,2)=sr(k) cpp(k,icomp,3)=si(k) enddo call polz c generate plot-parameter file write(plop(8:8),'(i1)') icomp open(2,file=plop) write(2,*) '0, 4, 24, 15,', ns1, ',', ns2,', 0.8' write(2,8) einfile write(2,8) ausf(icomp) write(2,8) synt(icomp) write(2,8) rest(icomp) 8 format(a,a) close(2) write(6,*) 'a plot-parameter file ',plop,' was generated' write(4,*) 'a plot-parameter file ',plop,' was generated' write(6,*) write(4,*) 225 continue c ********************************************************************** c at this point the results for the U V W components are c stored in the arrays csp and cpp. c csp(i,icomp) contains the system parameters, syspar. c cpp contains explicit parameters as well as complex poles and residues. c cpp(i,icomp,1) are the explicit parameters like amplitude, delay, corner periods etc. c cpp(i,icomp,2) are the real parts of complex poles c cpp(i,icomp,2) are the imag parts of complex poles c now transform parameters back to X, Y, Z in turn, and print the transformed values. call backtr(csp,cpp) do 226 icomp=1,3 write(6,104) write(4,104) write(6,104) write(4,104) write(6,110) xyz(icomp) write(4,110) xyz(icomp) 110 format(' final system parameters for the ',a1,' component') iter=0 q=zero write(6,103) (name(k),k=1,m) write(4,103) (name(k),k=1,m) write(6,104) iter,dsqrt(q),(cpp(k,icomp,1),k=1,m) write(4,104) iter,dsqrt(q),(cpp(k,icomp,1),k=1,m) syspar(0)=csp(0,icomp) do k=1,msys syspar(k)=csp(k,icomp) sr(k)=cpp(k,icomp,2) si(k)=cpp(k,icomp,3) enddo call polz 226 continue write(6,*) 'this protocol was saved as file trical.out' close(4) stop end subroutine inpar(m,maxit,qac,finac,init,dt,ns1,ns2) c read parameters from the parameter file parameter (mpar=12,msys=25) implicit double precision (a-h,o-z) character*55 zeile character*3 name,nam,typ,type logical nexttyp common /par/ x0(msys),rho(msys),x00(mpar),r00(mpar),typ(msys), & name(mpar) common /inp/ m0i,m1i,m2i,m10,m20 common /ali/ mali,malias,alias,mal data zero,one,two/0.d0,1.d0,2.d0/ ipar=0 nexttyp=.true. type=' ' open(3,file='trical.par',status='old') 1 format(a) c ****************************** read inversion parameters read(3,1) zeile do izeile=1,4 201 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 201 enddo 2 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 2 read(zeile,*) alias 3 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 3 read(zeile,*) m 4 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 4 read(zeile,*) m0i 5 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 5 read(zeile,*) m1i 6 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 6 read(zeile,*) m2i 7 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 7 read(zeile,*) maxit 8 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 8 read(zeile,*) qac 9 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 9 read(zeile,*) finac 10 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 10 read(zeile,*) ns1 11 read(3,1) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 11 read(zeile,*) ns2 if(m.eq.0) maxit=0 if(finac.gt.one) finac=one if(init.ge.0) then c ******************************************* read system parameters write(6,1) ' reading start parameters:' write(4,1) ' reading start parameters:' do 154 l=1,msys 23 read(3,1,end=14) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 23 if(zeile(1:3).eq.'END'.or.zeile(1:3).eq.'end') goto 14 if(l.gt.4.and.nexttyp) then read(zeile,'(a3)') type write(6,'(4x,a4)') type//':' write(4,'(4x,a4)') type//':' if(type.ne.'lp1'.and.type.ne.'hp1'.and. & type.ne.'lp2'.and.type.ne.'bp2'.and.type.ne.'hp2') then write(6,31) type write(4,31) type 31 format('subsystem ',a3,' not recognized - must be lp1, hp1, lp &2 ,bp2, or hp2') stop endif read(3,1,end=14) zeile if(zeile(1:1).eq.' '.or.zeile(1:1).eq.'') goto 23 if(zeile(1:3).eq.'END'.or.zeile(1:3).eq.'end') goto 14 endif read(zeile,*) nam,x0(l),rho(l) if(l.eq.1.and.nam.ne.'amp') then write(6,'(a)') 'first name must be amp' write(4,'(a)') 'first name must be amp' stop 1 else if(l.eq.2.and.nam.ne.'del') then write(6,'(a)') 'second name must be del' write(4,'(a)') 'second name must be del' stop 2 else if(l.eq.3.and.nam.ne.'sub') then write(6,'(a)') 'third name must be sub' write(4,'(a)') 'third name must be sub' stop 3 else if(l.eq.4.and.nam.ne.'til') then write(6,'(a)') 'fourth name must be til' write(4,'(a)') 'fourth name must be til' stop 4 endif write(6,16) nam,x0(l),rho(l),type write(4,16) nam,x0(l),rho(l),type 16 format(4x,a3,2x,2f10.3,3x,a3) if( (x0(2).ne.zero.or.rho(2).ne.zero).and. & (x0(3).ne.zero.or.rho(3).ne.zero) ) then write(6,17) write(4,17) 17 format(/" Parameters DEL and SUB should not be used at the same ti &me" & /" because a half-bridge does not delay the direct signal"/) c stop endif if(rho(l).ne.zero) then ipar=ipar+1 name(ipar)=nam x00(ipar)=x0(l) r00(ipar)=rho(l) endif if(type(3:3).eq.'2') nexttyp=.not.nexttyp typ(l)=type 154 continue write(6,101) write(4,101) 101 format(' parameterfile was not read to end') stop 14 if(ipar.eq.m.and.ipar.le.mpar.and.l.eq.5+m1i+2*m2i) goto 18 write(6,102) write(4,102) 102 format(' Parameter file is inconsistent. Check active parameters a &nd subsystems.') stop else c *************************** determine order of anti-alias filter write(6,15) write(6,13) m,m0i,m1i,m2i,maxit,qac,finac,alias 13 format(5i5,3f10.6) write(4,15) write(4,13) m,m0i,m1i,m2i,maxit,qac,finac,alias 15 format(' control parameters for the iteration:'/ &' m m0 m1 m2 maxit qac finac alias') m10=m1i m20=m2i m=0 m0i=0 m1i=0 m2i=0 do i=1,4 x0(i)=zero rho(i)=zero enddo x0(1)=one l=5 if(alias.lt.4.d0*dt) then write(6,103) alias write(4,103) alias 103 format(' specified anti-alias corner period is too small:',f10.3) stop else malias=int(6.d0/dlog10(alias/dt)+one) mali=malias/2 mal=malias-2*mali write(6,104) malias,alias write(4,104) malias,alias 104 format(' anti- alias- filter of order',i2,', T=',f6.2,' s') endif c ************************** anti-alias filter has order malias endif 18 if(mal.eq.1) then c anti-alias filter has odd order. Store first-order filter module as c parameter #5 after moving first-order system modules, if present, c to #6 and higher do 203 i=l-1,5,-1 x0(i+1)=x0(i) rho(i+1)=rho(i) typ(i+1)=typ(i) 203 continue l=l+1 m1i=m1i+1 x0(5)=alias rho(5)=zero typ(5)='lp1' endif wi=two*datan(one)/malias m2i=m2i+mali do 214 i=1,mali l2i=l+2*i-2 x0(l2i)=alias rho(l2i)=zero typ(l2i)='lp2' x0(l2i+1)=dsin(wi*(2*i-1)) rho(l2i+1)=zero typ(l2i+1)='lp2' 214 continue close(3) return end subroutine input(name,n,dt,x,tmin,tsec) dimension ix(60) double precision dt,x(n),tmin,tsec character iform*20,name*50,text*72,zeile*72 character zeil*120,zarr(120),b,code1*6,code2*6 logical seife b=' ' do j=1,n x(j)=0. enddo write(6,*) 'Opening file ',trim(name) write(4,*) 'Opening file ',trim(name) open(7,file=name,status='old') read(7,'(a)') zeile read(7,'(a)') zeile seife=index(zeile,'%').gt.0.or.index(zeile,'(').gt.0 if(seife) then write(6,*) 'file ',trim(name),' assumed to be in SEIFE format' else write(6,*) 'file ',trim(name),' assumed to be in ASL format' endif close(7) if(seife) then c read data in SEIFE format open(7,file=name,status='old') read(7,'(a)') text write(6,*) ' header: ',trim(text) 21 read(7,'(a)') zeile if(zeile(1:1).eq.'%') goto 21 20 read(zeile,1) nn,iform,dt,tmin,tsec 1 format(i10,a20,3f10.3) if(nn.gt.n) then write(*,*) 'sorry, too many data points. can handle only ', n write(*,*) 'for more, edit and recompile the source code' else n=nn endif write(6,22) n,trim(name) 22 format(' reading',i7,' samples from file ',a) if(iform(1:3).eq.'fre'.or.index(iform,'i').gt.0) then read(7,*,err=25,end=23) (x(j),j=1,n) else read(7,iform,err=25,end=23) (x(j),j=1,n) endif close(7) return else c read data in ASL format (such as written by Quanterra's Cimarron) open(7,file=name,status='old') read(7,'(a)') zeil text=zeil(1:72) write(6,*) 'header: ',trim(text) read(zeil,'(120a)') zarr nn=0 do i=1,119 if(zarr(i).ne.b.and.zarr(i+1).eq.b) then nn=nn+1 ix(nn)=i endif enddo read(zeil(1:ix(1)),'(a)') code1 read(zeil(ix(1)+2:ix(2)),'(a)') code2 read(zeil(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec read(zeil(ix(7)+2:ix(11)),*) t1000,srate,nn thr=ithr tmin=itmin tsec=itsec dt=1./srate tmin=tmin+60.*thr tsec=tsec+t1000/1000. if(nn.gt.n) then write(*,*) 'sorry, too many data points. can handle only ', n write(*,*) 'for more, edit and recompile the source code' else n=nn endif dt=1./srate tmin=tmin+60.*thr tsec=tsec+t100/100. c read(7,*,err=25,end=23) (x(j),j=1,n) c free-format read doesn't always work - compiler bug? c we must write our own routine nz=0 nval=0 do write(zeil,'(120x)') read(7,'(a)',end=31) zeil nz=nz+1 read(zeil,'(120a)') zarr nn=0 do i=1,119 if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1 enddo nn=min(nn,n-nval) read(zeil,*,err=32) (ix(j),j=1,nn) do j=1,nn x(nval+j)=ix(j) enddo nval=nval+nn enddo 31 write(6,*) nval,' samples read from file ',trim(name) write(4,*) nval,' samples read from file ',trim(name) close(7) return 32 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,nz,' of file ',trim(name) stop endif 23 n=j-1 write(6,24) n write(3,24) n 24 format('end of file after ',i8,' samples') close(7) return 25 write(6,26) j write(3,26) j 26 format(' Input error at sample # ',i8) stop end subroutine output(name,text,n,dt,x,tmin,tsec) double precision dt,x(n),xmax,tmin,tsec character name*5,text*(*),form*20 xmax=0.d0 do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) nvor=log10(max(xmax,2.d0))+1. ndec=max(0,10-nvor) write(form,3) ndec 3 format('(5f13.',i1,')') write(6,101) name write(4,101) name 101 format(' writing file ',a) open(8,file=name) write(8,'(a)') text write(8,1) n,form,dt,tmin,tsec 1 format(i10,a20,f10.6,2f10.3) write(8,form) (x(j),j=1,n) close(8) return end subroutine congrd(x,q,gnorm,iter,ein,aus,sum,einf,sta) c Method of conjugate gradients according to Fletcher und Reeves (1964) parameter (mpar=12,ns=480000) implicit double precision (a-h,o-z) dimension ein(ns),aus(ns),sum(ns),einf(ns),sta(ns) dimension x(mpar),g(mpar),d(mpar) dimension dd(mpar) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m data zero,two/0.d0,2.d0/ c partial derivatives do 1 k=1,m x(k)=x(k)+finac call quadr(x,qq,ein,aus,sum,einf,sta) x(k)=x(k)-finac-finac call quadr(x,qqq,ein,aus,sum,einf,sta) x(k)=x(k)+finac 1 g(k)=(qq-qqq)/two/finac c determine new direction of descent if(mod(iter-1,m).eq.0) then write(6,6) write(4,6) 6 format(1x) gnorm=0.d0 do 2 k=1,m gnorm=gnorm+g(k)**2 2 d(k)=-g(k) else ga=gnorm gnorm=zero do 3 k=1,m 3 gnorm=gnorm+g(k)**2 beta=gnorm/ga do 4 k=1,m 4 d(k)=-g(k)+beta*d(k) endif dlen=0.d0 do 5 k=1,m 5 dlen=dlen+d(k)**2 dlen=dsqrt(dlen) do 8 k=1,m 8 dd(k)=d(k)/dlen c search for minimum call mini(x,q,dd,ein,aus,sum,einf,sta) return end subroutine mini(x,q,d,ein,aus,sum,einf,sta) c new version, Oct. 2008, equivalent to the old version c x and q entered as start values, replaced by values of minimum c step size is between 1 and finac parameter (mpar=12,ns=480000) implicit double precision (a-h,o-z) logical exit3,exit4,loop4,interp dimension ein(ns),aus(ns),sum(ns),einf(ns),sta(ns) dimension x(mpar),xl(mpar),xm(mpar),xr(mpar),d(mpar) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m data one,two,eight /1.d0,2.d0,8.d0/ call move(xm,qm,x,q) do k=1,m xl(k)=x(k)-step*d(k) xr(k)=x(k)+step*d(k) enddo call quadr(xl,ql,ein,aus,sum,einf,sta) call quadr(xr,qr,ein,aus,sum,einf,sta) if(ql.lt.qm.and.qr.lt.qm) then c qm is the maximum of ql, qm, qr. This should never happen. stop 'Something went wrong. Try different search ranges.' endif loop_3: do if(qm.le.ql.and.qm.le.qr) then c qm is the minimum of ql, qm, qr. We can leave loop 3. if(step.lt.eight*finac) then exit3=.true. loop4=.false. interp=.true. else exit3=.false. step=step/eight do k=1,m xl(k)=xm(k)-step*d(k) xr(k)=xm(k)+step*d(k) enddo call quadr(xl,ql,ein,aus,sum,einf,sta) call quadr(xr,qr,ein,aus,sum,einf,sta) endif else exit3=.true. loop4=.true. if(ql.lt.qr) then c! reverse direction so we don't have to deal with two cases call move(x,q,xl,ql) call move(xl,ql,xr,qr) call move(xr,qr,x,q) do k=1,m d(k)=-d(k) enddo endif endif if(exit3) exit loop_3 enddo loop_3 if(loop4) then loop_4: do if(step.gt.one) then c! don't interpolate exit4=.true. interp=.false. call move(x,q,xr,qr) axi=step else c! double interval beyond xr step=two*step call move(xm,qm,xr,qr) do k=1,m xr(k)=xm(k)+step*d(k) enddo call quadr(xr,qr,ein,aus,sum,einf,sta) if(qr.lt.qm) then exit4=.false. else exit4=.true. interp=.true. endif endif if(exit4) exit loop_4 enddo loop_4 endif if(interp) then xi=(ql-qr)/(ql-two*qm+qr)/two*step axi=dabs(xi) do k=1,m x(k)=xm(k)+xi*d(k) enddo call quadr(x,q,ein,aus,sum,einf,sta) endif return end subroutine mini_s(x,q,d,ein,aus,sum,einf,sta) c original 'spaghetti' version of subroutine mini (of around 1985) c x and q entered as start values, replaced by values of minimum c step size is between 1 and finac parameter (mpar=12,ns=480000) implicit double precision (a-h,o-z) dimension ein(ns),aus(ns),sum(ns),einf(ns),sta(ns) dimension x(mpar),xl(mpar),xm(mpar),xr(mpar),d(mpar) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m data one,two,eight /1.d0,2.d0,8.d0/ call move(xm,qm,x,q) do 1 k=1,m xl(k)=x(k)-step*d(k) 1 xr(k)=x(k)+step*d(k) call quadr(xl,ql,ein,aus,sum,einf,sta) call quadr(xr,qr,ein,aus,sum,einf,sta) c! maximum encountered 3 if(ql.lt.qm.and.qr.lt.qm) then stop 222 endif c! go on if(ql.lt.qm.or.qr.lt.qm) goto 2 c! interval too small 11 if(step.lt.eight*finac) goto 8 c! divide interval step=step/eight do 7 k=1,m xl(k)=xm(k)-step*d(k) 7 xr(k)=xm(k)+step*d(k) call quadr(xl,ql,ein,aus,sum,einf,sta) call quadr(xr,qr,ein,aus,sum,einf,sta) goto 3 2 if(ql.lt.qr) then c! reverse direction call move(x,q,xl,ql) call move(xl,ql,xr,qr) call move(xr,qr,x,q) do 9 k=1,m 9 d(k)=-d(k) endif 4 if(step.gt.one) then c! don't interpolate call move(x,q,xr,qr) axi=step return else step=two*step call move(xm,qm,xr,qr) c! double interval do 5 k=1,m 5 xr(k)=xm(k)+step*d(k) call quadr(xr,qr,ein,aus,sum,einf,sta) if(qr.lt.qm) goto 4 c! again? endif 8 xi=(ql-qr)/(ql-two*qm+qr)/two*step axi=dabs(xi) do 6 k=1,m 6 x(k)=xm(k)+xi*d(k) call quadr(x,q,ein,aus,sum,einf,sta) return end subroutine move(x2,q2,x1,q1) parameter (mpar=12) implicit double precision (a-h,o-z) dimension x1(mpar),x2(mpar) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m do 1 k=1,m 1 x2(k)=x1(k) q2=q1 return end subroutine quadr(x,quad,ein,aus,sum,einf,sta) c This routine calculates the synthetic output signal an subtracts c it from the observed one, in order to determine the energy of c the residual which will be minimized with the conjugate-gradient c method. The energy is normalized to the energy of the signal. c The system is modelled as a parallel arrangement (sum) of first-order c and second-order filters, not as an arrangement in series (product). c For details see the publication by Schuessler. parameter (mpar=12,ns=480000,msys=25) implicit double precision (a-h,o-z) character name*3,typ*3 dimension x(mpar) dimension ein(ns),aus(ns),sum(ns),einf(ns),sta(ns) common /aux/ step,finac,axi,qn,dt,ns1,ns2,mq,init,m common /sys/ syspar(0:msys),sr(msys),si(msys),resid(msys), & residi(msys),m0,m1,m2 common /par/ x0(msys),rho(msys),x00(mpar),r00(mpar),typ(msys), & name(mpar) data zero,half,grav/0.d0,5.d-1,9.81d-3/ mq=mq+1 c write(*,*) 'quadr Aufruf ', mq ns21=ns2-ns1+1 c in order to be able to delay and advance the synthetics, we must c compute a few extra points outside the interval (ns1,ns2) ns14=ns1-24 ns24=ns2+24 call sysdef(x,dt,delay,compens1,compens2) if(dabs(delay/dt).ge.24d0) then write(6,101) write(4,101) 101 format(' Delay exceeds 24*dt.'/ & ' Use a smaller search range for the delay parameter.'/ & ' A range smaller than the sampling interval is recommended.'/ & ' If this does not help, then you must edit the time series.'/ & ' Remove an appropriate number of samples from the beginning of'/ & ' either the input or the output signal to reduce the delay.') stop endif call partl do 455 i=1,ns24 sum(i)=zero 455 continue ishift=ceiling(delay/dt) tau=delay-dt*ishift do 462 j=1,m1 i=m0+j call rekf1(dt*sr(i),resid(i)*dexp(-tau*sr(i)),ns24,ein,sum) 462 continue do 467 j=1,m2 i=m0+m1+j if(si(i).eq.zero) then call rekf1(dt*sr(i),resid(i)*dexp(-tau*sr(i)),ns24,ein,sum) call rekf1(dt*sr(i+m2),resid(i+m2)*dexp(-tau*sr(i+m2)),ns24,ein, & sum) else sre=-tau*sr(i) sim=-tau*si(i) zabs=dexp(sre) zre=zabs*dcos(sim) zim=zabs*dsin(sim) rre=resid(i)*zre-residi(i)*zim rim=resid(i)*zim+residi(i)*zre call rekf2(dt*sr(i),dt*si(i),rre,rim,ns24,ein,sum) endif 467 continue if(init.lt.0) return if(ishift.gt.0) then do 482 i=ns2,ns1,-1 sum(i)=sum(i-ishift) 482 continue else if(ishift.lt.0) then do 483 i=ns1,ns2 sum(i)=sum(i-ishift) 483 continue endif c sum is now the synthetics. Correct for galvanic coupling or c shake-table tilt if required, form difference with observed c output, and remove offset. c sub is the fraction of the input signal to be subtracted c from the output. c til ist the tilt of the shake table in microradians per c millimeter of motion. grav is the gravity in km/s^2. c Note: motion and tilt as seen by the sensor, so if the c motion is transverse, the tilt parameter may come out c very large! if(compens2.ne.zero) then sta(1)=half*sum(ns1) do i=ns1+1,ns2 ii=i-ns1+1 sta(ii)=sta(ii-1)+half*(sum(i-1)+sum(i)) enddo call polytrend(1,sta,ns21) suma=sta(1) sta(1)=half*suma do i=ns1+1,ns2 ii=i-ns1+1 sumn=sta(ii) sta(ii)=sta(ii-1)+half*(suma+sumn) suma=sumn enddo call polytrend(2,sta,ns21) endif comp=compens2*grav*dt**2 do i=ns1,ns2 ii=i-ns1+1 sum(i)=aus(i)-sum(i)-compens1*einf(i)+comp*sta(ii) enddo do i=ns1,ns2 ii=i-ns1+1 sta(ii)=sum(i) enddo call polytrend(3,sta,ns21) c call polytrend(1,sta,ns21) do i=ns1,ns2 ii=i-ns1+1 sum(i)=sta(ii) enddo offset=zero qua=zero do i=ns1,ns2 offset=offset+sum(i) enddo offset=offset/ns21 do i=ns1,ns2 sum(i)=sum(i)-offset qua=qua+sum(i)**2 enddo quad=qua/qn c At the end of the program, sum is the misfit (with its c mean removed), not the synthetic. quad is the normalized c energy of the misfit, as a fraction of the total energy. return end subroutine sysdef(x,dt,delay,compens1,compens2) parameter (mpar=12,msys=25) implicit double precision (a-h,o-z) dimension x(mpar),xx(msys) character*3 typ,name common /par/ x0(msys),rho(msys),x00(mpar),r00(mpar),typ(msys), & name(mpar) common /sys/ syspar(0:msys),sr(msys),si(msys),resid(msys), & residi(msys),m0,m1,m2 common /inp/ m0i,m1i,m2i,m10,m20 data zero,one /0.d0,1.d0/ C Definition of the transfer function (Laplace Transform): C syspar(0 : m0): coefficients of the nominator polynomial, order m0 C syspar(m0+1 : m0+m1) corner freq. of first-order partials C syspar(m0+m1+1 : m0+m1+m2) corner freq. of second-order partials C syspar(m0+m1+m2+1 : m0+m1+2*m2) numerical damping factors C C C Sum syspar(i) * s**i C i=0,m0 C T(s) = -------------------------------------------------------- C Product (s+w1(i))*Product(s**2+2*s*w2(i)*h(i)+w2(i)**2) C i=1,m1 i=1,m2 C C C where: w1(i) = syspar(m0+i) , i = 1 ... m1 C w2(i) = syspar(m0+m1+i) , i = 1 ... m2 C h(i) = syspar(m0+m1+m2+i) , i = 1 ... m2 C In addition the signal may be delayed by 'delay' (abs < +-dt), C diminished by a fraction 'compens1' of the input signal to compensate C for the voltage coupled into the coil of geophones C by the calibration current, and a fraction 'compens2' of the twice- C integrated signal may be removed in order to compensate for the tilt of C shake tables. C current system parameters (order is different from parameter file!) zpi=8.d0*datan(one) m0=m0i m1=m1i m2=m2i mp=4+m1+2*m2 ipar=0 do 546 i=1,mp if(rho(i).gt.zero) then ipar=ipar+1 xx(i)=x0(i)+x(ipar)*rho(i) c if (icount.lt.3) write(6,547) i,name(ipar),x0(i),xx(i),rho(i) else xx(i)=x0(i) c if (icount.lt.3) write(6,547) i,'***',x0(i),xx(i),rho(i) endif 547 format(i2,2x,a3,3f10.3) 546 continue c icount=icount+1 c if(icount.ge.3) stop gain=xx(1)*dt delay=xx(2) compens1=xx(3) compens2=xx(4) c Numerator factor for first-order subsystems do 563 i=5,m1+4 c write(6,*) i,typ(i) if(typ(i).eq.'hp1') then m0=m0+1 else if(typ(i).eq.'lp1') then gain=gain*zpi/xx(i) else write(6,102) typ(i) write(4,102) typ(i) 102 format('wrong type: ',a) stop endif 563 continue c Numerator factor for second-order subsystems do 571 i=m1+5,m1+2*m2+3,2 c write(6,*) i,typ(i) if(typ(i).eq.'hp2') then m0=m0+2 else if (typ(i).eq.'bp2') then m0=m0+1 gain=gain*zpi/xx(i) else if (typ(i).eq.'lp2') then gain=gain*(zpi/xx(i))**2 else write(6,102) typ(i) write(4,102) typ(i) stop endif 571 continue if(m0.gt.m1+2*m2-2) then write(6,105) write(4,105) 105 format(' m0 is too large - stop') stop endif if(m0.lt.-1) then write(6,103) write(4,103) 103 format(' m0 is too small - stop') stop endif if(m0.ge.0) then do 594 i=0,m0-1 syspar(i)=zero 594 continue syspar(m0)=gain do 598 i=1,m1 syspar(m0+i)=zpi/xx(i+4) 598 continue else m0=0 syspar(m0)=gain do 604 i=1,m1 syspar(i)=zpi/xx(i+4) 604 continue m1=m1+1 syspar(m1)=zero endif do 611 i=1,m2 syspar(m0+m1+i)=zpi/xx(2*i+m1+3) syspar(m0+m1+m2+i)=xx(2*i+m1+4) 611 continue return end subroutine polz parameter (msys=25) implicit double precision (a-h,o-z) common /inp/ m0i,m1i,m2i,m10,m20 common /sys/ syspar(0:msys),sr(msys),si(msys),resid(msys), & residi(msys),m0,m1,m2 common /ali/ mali,malias,alias,mal data zero,one/0.d0,1.d0/ C first-order partials: poles sr and residues resid are real and have c indices m0+i with i=1 .. m1 C second-order partials: poles (sr,si) and residues (resid,residi) are c pairwise complex-conjugate and stored in c sr(m0+m1+i), si(m0+m1+i) and sr(m0+m1+m2+i), si(m0+m1+m2+i), i=1 .. m2 write(6,*) write(4,*) write(6,*) 'Transfer function (counts in ==> counts out):' write(4,*) 'Transfer function (counts in ==> counts out):' write(6,*) write(4,*) if(m0.eq.1) then write(6,*) 'The Laplace transform has 1 zero at zero frequency' write(4,*) 'The Laplace transform has 1 zero at zero frequency' else write(6,*) 'The Laplace transform has',m0, & ' zeros at zero frequency' write(4,*) 'The Laplace transform has',m0, & ' zeros at zero frequency' endif write(6,*) write(4,*) if(m10.eq.1) then write(6,*) 'It has 1 real pole, in rad/s:' write(4,*) 'It has 1 real pole, in rad/s:' else write(6,*) 'It has',m10,' real poles, in rad/s:' write(4,*) 'It has',m10,' real poles, in rad/s:' endif c mal=malias-2*mali do i=1,m10 write(6,*) sngl(sr(m0+i+mal)) write(4,*) sngl(sr(m0+i+mal)) enddo write(6,*) write(4,*) mm=m20 if(mm.eq.1) then write(6,*) 'and 1 pair of complex-conjugate poles, in rad/s:' write(4,*) 'and 1 pair of complex-conjugate poles, in rad/s:' else write(6,*) 'and',mm,' pairs of complex-conjugate poles, in rad/s:' write(4,*) 'and',mm,' pairs of complex-conjugate poles, in rad/s:' endif do i=1,m20 write(6,*) sngl(sr(m0+m1+i)), sngl(si(m0+m1+i)) write(6,*) sngl(sr(m0+m1+i)), -sngl(si(m0+m1+i)) write(4,*) sngl(sr(m0+m1+i)), sngl(si(m0+m1+i)) write(4,*) sngl(sr(m0+m1+i)), -sngl(si(m0+m1+i)) enddo write(6,*) write(4,*) return end subroutine partl parameter (msys=25) implicit double precision (a-h,o-z) common /sys/ syspar(0:msys),sr(msys),si(msys),resid(msys), & residi(msys),m0,m1,m2 data zero,one/0.d0,1.d0/ C factorization of the transfer function. C poles und residues (real and imag.parts) stored like system parameters C first-order partials: poles sr and residues resid are real and have c indices m0+i with i=1 .. m1 C second-order partials: poles (sr,si) and residues (resid,residi) are c pairwise complex-conjugate and stored in c sr(m0+m1+i), si(m0+m1+i) and sr(m0+m1+m2+i), si(m0+m1+m2+i), i=1 .. m2 C residues are stored in resid(m0+m1+i), residi(mo+m1+i), i=1 .. 2*m2 C determine the complex poles do 639 i=1,m1 sr(m0+i)=-syspar(m0+i) si(m0+i)=zero 639 continue do 643 i=1,m2 w=syspar(m0+m1+i) h=syspar(m0+m1+m2+i) if(dabs(h-one).lt.1.d-6) then write(6,111) i write(4,111) i 111 format(/' SORRY, the program fails when damping is exactly', & ' critical (=1.000000)'/' Try again with a different start ', & 'value for damping #',i2) stop endif if(h.lt.zero) then write(6,101) i write(4,101) i 101 format(' Negative damping encountered. Reduce search range for & damping #',i2) stop endif wh=w*h if(h.lt.one) then whh=w*dsqrt(one-h**2) sr(m0+m1+i)=-wh si(m0+m1+i)=whh sr(m0+m1+m2+i)=-wh si(m0+m1+m2+i)=-whh else whh=w*dsqrt(h**2-one) sr(m0+m1+i)=-wh+whh si(m0+m1+i)=zero sr(m0+m1+m2+i)=-wh-whh si(m0+m1+m2+i)=zero endif 643 continue C calculate residues do 1 i=m0+1,m0+m1+2*m2 cresr=syspar(0) cresi=zero spowr=one spowi=zero do 2 k=1,m0 spowre=spowr*sr(i)-spowi*si(i) spowi=spowr*si(i)+spowi*sr(i) spowr=spowre cresr=cresr+syspar(k)*spowr 2 cresi=cresi+syspar(k)*spowi do 3 k=m0+1,m0+m1+2*m2 if(k.ne.i) then c cres=cres/(s(i)-s(k)) sdr=sr(i)-sr(k) sdi=si(i)-si(k) sdq=sdr*sdr+sdi*sdi cresre=(cresr*sdr+cresi*sdi)/sdq cresi=(cresi*sdr-cresr*sdi)/sdq cresr=cresre endif 3 continue resid(i)=cresr 1 residi(i)=cresi return end subroutine rekf1(dt,r,n,ein,sum) parameter(ns=480000) implicit double precision (a-h,o-z) dimension ein(ns),sum(ns) z=dexp(dt) aus0=r*ein(1) sum(1)=sum(1)+aus0 do 691 i=2,n aus0=r*ein(i)+z*aus0 sum(i)=sum(i)+aus0 691 continue return end subroutine rekf2(sre,sim,rre,rim,n,ein,sum) parameter(ns=480000) implicit double precision (a-h,o-z) dimension ein(ns),sum(ns) data two/2.d0/ zabs=dexp(sre) zre=zabs*dcos(sim) zim=zabs*dsin(sim) f0=two*rre f1=-two*(rre*zre+rim*zim) g1=two*zre g2=-zabs**2 aus0=f0*ein(1) sum(1)=sum(1)+aus0 aus1=aus0 aus0=f0*ein(2)+f1*ein(1)+g1*aus1 sum(2)=sum(2)+aus0 do 719 i=3,n aus2=aus1 aus1=aus0 aus0=f0*ein(i)+f1*ein(i-1)+g1*aus1+g2*aus2 sum(i)=sum(i)+aus0 719 continue return end subroutine polytrend(mm,x,n) c remove polynomial trend parameter(ndi=8,ns=480000) implicit real*8 (a-h,o-z) dimension b(ndi),c(ndi,ndi),a(ndi),x(ns) m=min(mm,ndi-1) fnh=n/2. zero=0.d0 one=1.d0 do j=1, m+1 do k=1, m+1 c(j,k)=zero do i=1, n c(j,k)=c(j,k)+(dble(i)/fnh-one)**(j+k-2) enddo enddo b(j)=zero do i=1,n b(j)=b(j)+(dble(i)/fnh-one)**(j-1)*x(i) enddo enddo call gauss(c,m+1,ndi,b,a) 100 format(i5,e15.6) do i=1,n xpol=a(m+1) do j=m,1,-1 xpol=xpol*(dble(i)/fnh-one)+a(j) enddo x(i)=x(i)-xpol enddo return end subroutine gauss(aik,m,n,rs,f) c solve linear equations implicit real*8 (a-h,o-z) parameter(ndi=8,ndi1=ndi+1) dimension aik(n,n),rs(n),f(n),h(ndi1),imax(ndi) do 1401 j=1,m aikmax=0.d0 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 subroutine trans(aus3,n) parameter(ns=480000) implicit real*8 (a-h,o-z) dimension aus3(ns,3) w2=dsqrt(2.d0) w3=dsqrt(3.d0) w6=dsqrt(6.d0) aux=-2.d0/w6 auy=0.d0 auz=w2/w6 avx=1.d0/w6 avy=w3/w6 avz=auz awx=avx awy=-avy awz=auz do i=1,n u=aux*aus3(i,1)+auy*aus3(i,2)+auz*aus3(i,3) v=avx*aus3(i,1)+avy*aus3(i,2)+avz*aus3(i,3) w=awx*aus3(i,1)+awy*aus3(i,2)+awz*aus3(i,3) aus3(i,1)=u aus3(i,2)=v aus3(i,3)=w enddo write(4,*) '==> signals were transformed from X Y Z to U V W' write(6,*) '==> signals were transformed from X Y Z to U V W' return end subroutine backtr(csp,cpp) parameter (msys=25) implicit real*8 (a-h,o-z) dimension csp(0:msys,3),cpp(msys,3,3) data zero,eins,zwei,drei,sechs /0.d0,1.d0,2.d0,3.d0,6.d0/ axu=zwei/drei axv=eins/sechs axw=axv ayu=zero ayv=eins/zwei ayw=ayv azu=eins/drei azv=azu azw=azu do i=0,msys x=axu*csp(i,1)+axv*csp(i,2)+axw*csp(i,3) y=ayu*csp(i,1)+ayv*csp(i,2)+ayw*csp(i,3) z=azu*csp(i,1)+azv*csp(i,2)+azw*csp(i,3) csp(i,1)=x csp(i,2)=y csp(i,3)=z enddo do k=1,3 do i=1,msys x=axu*cpp(i,1,k)+axv*cpp(i,2,k)+axw*cpp(i,3,k) y=ayu*cpp(i,1,k)+ayv*cpp(i,2,k)+ayw*cpp(i,3,k) z=azu*cpp(i,1,k)+azv*cpp(i,2,k)+azw*cpp(i,3,k) cpp(i,1,k)=x cpp(i,2,k)=y cpp(i,3,k)=z enddo enddo write(4,*) '==> parameters were transformed from U V W to X Y Z' write(6,*) '==> parameters were transformed from U V W to X Y Z' return end