C C======================================================================= SUBROUTINE datiris(sta,chan,ifl) C======================================================================= c CHARACTER*4 stal character*6 sta,chan character*40 ligne,fich character*3 comp dimension iuni(5),ia(5),juni(4),ja(4),jk(8) COMPLEX POLES, ZEROS character*1 tt,tq COMMON/HEADIR/A0,POLES(30,4),ZEROS(30,4),coen(60,5),coed(60,5), *dtp(4),dtc(5),nsp,npp(4),nzz(4),nsc,ncn(5),ncd(5),kuni,tt(4),tq(5) C DATA IU/20/,nscm/5/,npzm/30/,ncfm/60/,npm/4/ C open(iu,file='/usr/ipg/barbara/source/cali.seed', * status='old',form='formatted',iostat=it) if(it.ne.0) then write(*,*)'fichier calibration non trouve' nsp=0 nsc=0 return endif do 10 k=1,8 10 jk(k)=0 read(iu,'(a40)',end=800) ligne C 100 READ (ligne,2000) STAL,comp,A0,nsp,nsc if(nsp.gt.npm) go to 850 if(nsc.gt.nscm) go to 850 if(stal.ne.sta(1:4).or.comp.ne.chan(1:3)) then 110 read(iu,'(a)',end=800) ligne if(ligne(1:1).eq.' ') go to 110 go to 100 endif print*,'stal..',stal,comp,a0,nsp,nsc if (nsp.eq.0)go to 121 do 120 k=1,nsp READ (IU,'(2i5,2x,a1,2x,2i5,f5.2)') juni(k),ja(k),tt(k), * npp(k),nzz(k),dtp(k) print*,juni(k),ja(k),tt(k),npp(k),nzz(k),dtp(k) if(npp(k).gt.npzm) go to 850 if(nzz(k).gt.npzm) go to 850 if(tt(k).eq.'D'.and.dtp(k).lt..0001) go to 860 read(iu,2001) (POLES(I,k),I=1,npp(k)) READ (IU,2001) (ZEROS(I,k),I=1,nzz(k)) 120 continue 121 continue 2000 FORMAT (A4,2x,a3,E15.7,2i5) 2001 FORMAT ((4E17.5)) 2002 FORMAT ((6E17.5)) if(nsc.eq.0)go to 131 do 130 k=1,nsc READ (IU,'(2i5,2x,a1,2x,2i5,f5.2)') iuni(k),ia(k),tq(k), * ncn(k),ncd(k),dtc(k) print*,iuni(k),ia(k),tq(k),ncn(k),ncd(k),dtc(k) if(ncn(k).gt.ncfm) go to 850 if(ncd(k).gt.ncfm) go to 850 c if(tq(k).eq.'D'.and.dtc(k).lt..0001) go to 860 if(ncn(k).ne.0) then read(iu,2002) (coen(i,k),i=1,ncn(k)) endif if(ncd(k).ne.0) then read(iu,2002) (coed(i,k),i=1,ncd(k)) endif 130 continue 131 continue c c tests de compatibilite c poles et zeros a partie reelle negative c on inverse le signe completement c il faudrait tester sur l'existence du complexe conjugue do 150 k=1,nsp do 140 i=1,npp(k) if(real(poles(i,k)).gt.0.) then write(*,*)'correction pole: ',poles(i,k) poles(i,k)=-poles(i,k) endif 140 continue do 150 i=1,nzz(k) if(real(zeros(i,k)).gt.0.) then write(*,*)'correction zero: ',zeros(i,k) zeros(i,k)=-zeros(i,k) endif 150 continue c do 160 k=1,nsp ik=juni(k) jk(ik)=jk(ik)+1 ik=ja(k) 160 jk(ik)=jk(ik)+1 do 170 k=1,nsc ik=iuni(k) jk(ik)=jk(ik)+1 ik=ia(k) 170 jk(ik)=jk(ik)+1 c write(*,*)'test unite: ',(jk(k),k=1,8) do 180 k=1,8 if(mod(jk(k),2).ne.0) luni=luni+1 180 continue if(luni.ne.2) write(*,*) 'attention aux unites de sortie' if(mod(jk(4),2).eq.1) kuni=1 if(mod(jk(5),2).eq.1) kuni=2 if(mod(jk(6),2).eq.1) kuni=3 c c write(*,*) 'les parametres trouves pour ',sta,' sont:' c write(*,*) 'a0 = ',a0 c write(*,*) nsp, ' sequences de poles et zeros' c do 200 k=1,nsp c write(*,*) npp(k), 'poles', nzz(k),' zeros' c write(*,*) (poles(i,k),i=1,npp(k)) c write(*,*) (zeros(i,k),i=1,nzz(k)) c 200 continue c write(*,*) ' kuni=',kuni c write(*,*) nsc, ' sequences de coefficients' c write(*,*) (ncn(i),ncd(i),i=1,nsc) c close(iu) return c 800 WRITE(*,*)'PARAMETRES DE LA STATION NON TROUVES' go to 900 860 write(*,*)'attention a donner le pas pour les fitres digitaux' go to 900 850 write(*,*)'attention aux dimensions' 900 nsp =0 nsc=0 close(iu) RETURN END C C======================================================================= SUBROUTINE RIRIS (F,R,ifl,sta) C======================================================================= c COMPLEX s,r,anum COMPLEX POLES, ZEROS character*1 tt,tq character*6 sta COMMON/HEADIR/A0,POLES(30,4),ZEROS(30,4),coen(60,5),coed(60,5), *dtp(4),dtc(5),nsp,npp(4),nzz(4),nsc,ncn(5),ncd(5),kuni,tt(4),tq(5) C DATA PI2/6.2831852/ R=CMPLX(1.,0.) IF ((nsp.EQ.0).AND.(nsc.EQ.0)) RETURN c c reponse par poles et zeros c if(nsp.ne.0) then do 50 k=1,nsp c print*,'nsp k',nsp,k if(tt(k).eq.'A') then S=CMPLX(0.,F*PI2) else if(tt(k).eq.'B') then s=cmplx(0.,f) else if(tt(k).eq.'D') then c write(*,*) 'test',tt(k) aa=pi2*f*dtp(k) s=cmplx(cos(aa),sin(aa)) r=r*s**(npp(k)-nzz(k)) else write(*,*)'cas non defini' return endif endif endif IF (NPP(k).NE.0) THEN DO 20 I=1,NPP(k) 20 R=R/(S-POLES(I,k)) endif IF (NZZ(k).NE.0) THEN DO 10 I=1,NZZ(k) 10 R=R*(S-ZEROS(I,k)) endif 50 continue endif c c ajout des polynomes s'il faut C if(nsc.ne.0) then c print*,'nsc',nsc do 100 k=1,nsc c print*,'tq(k)ncn(k) ncd',tq(k),ncn(k),ncd(k) if(tq(k).ne.'D') then write(*,*)'attention cas de coefficients pour filtres' write(*,*)' analogiques non prevu' write(*,*)'le filtre',k,' n''est pas pris en compte' go to 100 else aa=pi2*f*dtc(k) s=cmplx(cos(aa),-sin(aa)) if(dtc(k).eq.0.)then c print*,k,dtc(k) s=cmplx(0.,f) endif if(ncn(k).ne.0) then anum=cmplx(0.,0.) do 110 j=1,ncn(k) 110 anum=anum+coen(j,k)*s**(j-1) r=r*anum endif c if(ncd(k).ne.0) then anum=cmplx(1.,0.) do 120 j=1,ncd(k) 120 anum=anum-coed(j,k)*s**(j) r=r/anum endif endif c 100 continue endif c print*,'a0',a0 R=A0*R c C sortie dans l'unite voulue C if(sta(1:3).eq.'NNA')kuni=2 if(sta(1:3).eq.'ESK')kuni=2 if(sta(1:3).eq.'PFO')kuni=2 print*,'kuni',kuni if(kuni.eq.0) return ip=ifl-kuni if(sta(1:3).eq.'NNA')ip=1 if(sta(1:3).eq.'ESK')ip=1 if(sta(1:3).eq.'PFO')ip=1 print* ,'ip',ip s=cmplx(0.,pi2*f) r=r/s**ip C RETURN END C