      SUBROUTINE APRFLT(FSPEC,AMPSQ,PHASE,NPOINTS,F1,F2)                        
C                                                                               
C  MODULE         APRFLT10                                                      
C  VERSION        1.0                                                           
C  LAST MODIFIED  20/1/82                                                       
C  AUTHOR         DAVE POWELL                                                   
C  ROUTINE NO     30                                                            
C                                                                               
C  THIS SUBROUTINE FINDS THE MAG(SQUARED)/PHASE RESPONSE OF A FILTER            
C  WHOSE SINGULARITIES ARE CONTAINED IN THE ARRAY 'FSPEC'. IF THE               
C  TYPE OF 'FSPEC' IS S-PLANE THEN F1 AND F2 ARE THE LOWER AND UPPER            
C  FREQUENCIES ,RESPECTIVELY , AT WHICH THE AMPLITUDE/PHASE RESPONSE            
C  IS EVALUATED . IF THE TYPE OF 'FSPEC' IS Z-PLANE THEN THE RESPONSE           
C  IS EVALUATED FROM F1*(SAMPLING FREQUENCY) TO F2*(SAMPLING FREQUENCY)         
C                                                                               
C  PARAMETERS                                                                   
C  **********                                                                   
C                                                                               
C  FSPEC    R*4A  ARRAY CONTAINING S OR Z PLANE SINGULARITIES                   
C  AMPSQ    R*4A  ARRAY TO BE FILLED WITH AMPLITUDE SQUARED RESPONSE            
C  PHASE    R*4A  ARRAY TO BE FILLED WITH PHASE RESPONSE                        
C  NPOINTS  I*4   NO OF POINTS AT WHICH MAG/PHASE TO BE EVALUATED               
C  F1       R*4   LOWEST FREQ AT WHICH MAG/PHASE TO BE EVALUATED                
C  F2       R*4   HIGHEST FREQ AT WHICH MAG/PHASE TO BE EVALUATED               
C                 IF TYPE OF 'FSPEC' IS Z-PLANE THEN RESPONSE IS                
C                 EVALUATED FROM F1*(SAMPLING FREQ) TO F2*(SAMPLING FREQ)       
C                                                                               
      INTEGER*4 NPOINTS                                                         
      REAL*4 FSPEC(9999),AMPSQ(NPOINTS),PHASE(NPOINTS),F1,F2                    
C                                                                               
C  EXTERNAL VARIABLES                                                           
C  ******************                                                           
C                                                                               
$INCLUDE SPLANE10.CMN(NLIST)                                                    
$INCLUDE ZPLANE10.CMN(NLIST)                                                    
C                                                                               
C  LOCAL VARIABLES                                                              
C  ***************                                                              
C                                                                               
C  IF       I*4   POINTER AT FREQ N0 (S & Z)                                    
C  WFREQ    R*4   ANG FREQ AT WHICH MAG/PHASE EVALUATED (S & Z)                 
C  IPS      I*4   POINTER AT SINGULARITY (S)                                    
C  IPZ      I*4   POINTER AT SINGULARITY (Z)                                    
C  SPOS     C*8   POSITION OF S-PLANE SINGULARITY (S)                           
C  ZPOS     C*8   POSITION OF Z-PLANE SINGULARITY (Z)                           
C  SGAIN    C*8   GAIN OF S-PLANE SINGULARITY (S)                               
C  ZGAIN    C*8   GAIN OF Z-PLANE SINGULARITY (Z)                               
C  LGAIN    R*4   MAG GAIN OF S-PLANE SINGULARITY(dB's) (S & Z)                 
C  PHGAIN   R*4   PHASE OF SINGULARITY (S & Z)                                  
C  TYPE     I*4   TYPE OF SINGULARITY (S & Z)                                   
C  PI       R*4   3.14159265....... (S & Z)                                     
C  INCFRQ   R*4   INCREMENTAL FREQUENCY (S & Z)                                 
C  CON      R*4   180.0/3.14159... (S & Z)                                      
C  DSQ      R*4   DISTANCE(SQUARED) FROM FREQ TO SINGULARITY (S & Z)            
C  RTN      I*4   ROUTINE NUMBER (S & Z)                                        
C  IER      I*4   ERROR NUMBER (S & Z)                                          
C  Z        C*8   POINT IN Z-PLANE AT WHICH RESPONSE EVALUATED (Z)              
C                                                                               
      REAL*4 WFREQ,LGAIN,CON,PI,INCFRQ,PHGAIN,DSQ                               
      INTEGER*4 IF,IPS,TYPE,IER,RTN,IPZ                                         
      COMPLEX*8 SPOS,SGAIN,ZPOS,ZGAIN,Z                                         
      PARAMETER(RTN=30)                                                         
C                                                                               
C                                                                               
C     (1) CHECK PARAMETERS                                                      
C     ********************                                                      
C                                                                               
C     Check type of FSPEC is either s or z plane.                               
      IF(INT(FSPEC(PSX)).NE.PSV.AND.INT(FSPEC(PZX)).NE.PZV)GOTO 9001            
C     Check NPOINTS is greater than one.                                        
      IF(NPOINTS.LE.1) GOTO 9002                                                
C     Check F1<F2                                                               
      IF (F1.GE.F2) GOTO 9003                                                   
C     Jump to Z-plane loops if 'FSPEC' type is Z-plane.                         
      IF(INT(FSPEC(PZX)).EQ.PZV) GOTO 4000                                      
C                                                                               
C     (2) EVALUATE VARIABLES (S-PLANE)                                          
C     ********************************                                          
C                                                                               
C                                                                               
C     Jump out if number of singularities less than or equal to 0               
      IF (INT(FSPEC(PSN)).LE.0) GOTO 9005                                       
      PI=4.0*ATAN(1.0)                                                          
      W1=2.0*PI*F1                                                              
      W2=2.0*PI*F2                                                              
      INCFRQ=(W2-W1)/FLOAT(NPOINTS-1)                                           
      CON=180.0/PI                                                              
C                                                                               
C     (3) FOR EACH FREQUENCY LOOP OVER SINGULARITIES (S-PLANE)                  
C     ********************************************************                  
C                                                                               
C     Loop over frequencies.                                                    
      DO 3001 IF=1,NPOINTS                                                      
        WFREQ=INCFRQ*FLOAT(IF-1)+W1                                             
        AMPSQ(IF)=0.0                                                           
        PHASE(IF)=0.0                                                           
C                                                                               
C       Loop over singularities.                                                
        DO 3000 IS=1,INT(FSPEC(PSN))                                            
          IPS=PSH+1+PSE*(IS-1)                                                  
          TYPE=INT(FSPEC(IPS+PST))                                              
          SPOS=CMPLX(FSPEC(IPS+PSR),FSPEC(IPS+PSI))                             
          SGAIN=CMPLX(FSPEC(IPS+PSG),FSPEC(IPS+PSJ))                            
C                                                                               
C         Real zero.                                                            
          IF (TYPE.EQ.PSTRZ) THEN                                               
            DSQ=REAL(SPOS)**2+WFREQ**2                                          
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            PHASE(IF)=PHASE(IF)+CON*ATAN2(WFREQ,-REAL(SPOS))                    
C                                                                               
C         Real pole.                                                            
          ELSEIF (TYPE.EQ.PSTRP) THEN                                           
            DSQ=REAL(SPOS)**2+WFREQ**2                                          
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            PHASE(IF)=PHASE(IF)-CON*ATAN2(WFREQ,-REAL(SPOS))                    
C                                                                               
C         Complex zero.                                                         
          ELSEIF (TYPE.EQ.PSTCZ) THEN                                           
            DSQ=REAL(SPOS)**2+(WFREQ-AIMAG(SPOS))**2                            
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(WFREQ-AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)+TEMP                                            
C                                                                               
C         Complex pole.                                                         
          ELSEIF (TYPE.EQ.PSTCP) THEN                                           
            DSQ=REAL(SPOS)**2+(WFREQ-AIMAG(SPOS))**2                            
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(WFREQ-AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)-TEMP                                            
C                                                                               
C         Complex conjugate zero pair.                                          
          ELSEIF (TYPE.EQ.PST2Z) THEN                                           
            DSQ=REAL(SPOS)**2+(WFREQ-AIMAG(SPOS))**2                            
            DSQ=DSQ*(REAL(SPOS)**2+(WFREQ+AIMAG(SPOS))**2)                      
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(WFREQ-AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)+TEMP                                            
            TEMP=CON*ATAN2(WFREQ+AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)+TEMP                                            
C                                                                               
C         Complex conjugate pole pair.                                          
          ELSEIF (TYPE.EQ.PST2P) THEN                                           
            DSQ=REAL(SPOS)**2+(WFREQ-AIMAG(SPOS))**2                            
            DSQ=DSQ*(REAL(SPOS)**2+(WFREQ+AIMAG(SPOS))**2)                      
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(WFREQ-AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)-TEMP                                            
            TEMP=CON*ATAN2(WFREQ+AIMAG(SPOS),-REAL(SPOS))                       
            PHASE(IF)=PHASE(IF)-TEMP                                            
C                                                                               
C         Type not known.                                                       
          ELSE                                                                  
            GOTO 9004                                                           
C                                                                               
          ENDIF                                                                 
C                                                                               
C         Add on magnitude and phase introduced by gain of singularity          
          LGAIN=20.0*ALOG10(CABS(SGAIN))                                        
          PHGAIN=CON*ATAN2(AIMAG(SGAIN),REAL(SGAIN))                            
          PHASE(IF)=PHASE(IF)+PHGAIN                                            
          AMPSQ(IF)=AMPSQ(IF)+LGAIN                                             
C                                                                               
3000    CONTINUE                                                                
3001  CONTINUE                                                                  
      RETURN                                                                    
C                                                                               
4000  CONTINUE                                                                  
C                                                                               
C     (4) EVALUATE VARIABLES (Z-PLANE)                                          
C     ********************************                                          
C                                                                               
C     Jump out if number of singularities less than or equal to 0               
      IF (INT(FSPEC(PZN)).LE.0) GOTO 9005                                       
      PI=4.0*ATAN(1.0)                                                          
      INCFRQ=(F2-F1)/FLOAT(NPOINTS-1)                                           
      CON=180.0/PI                                                              
C                                                                               
C                                                                               
C     (5) FOR EACH FREQ LOOP OVER SINGULARITIES (Z-PLANE)                       
C     ***************************************************                       
C                                                                               
C     Loop over frequencies.                                                    
      DO 5100 IF=1,NPOINTS                                                      
        WFREQ=INCFRQ*FLOAT(IF-1)+F1                                             
        AMPSQ(IF)=0.0                                                           
        PHASE(IF)=0.0                                                           
        Z=CEXP(CMPLX(0.0,WFREQ*PI/2.0))                                         
C                                                                               
C       Loop over singularities.                                                
        DO 5000 IZ=1,INT(FSPEC(PZN))                                            
          IPZ=PZH+1+PZE*(IZ-1)                                                  
          TYPE=INT(FSPEC(IPZ+PZT))                                              
          ZPOS=CMPLX(FSPEC(IPZ+PZR),FSPEC(IPZ+PZI))                             
          ZGAIN=CMPLX(FSPEC(IPZ+PZG),FSPEC(IPZ+PZJ))                            
C                                                                               
C         Real zero.                                                            
          IF (TYPE.EQ.PZTRZ) THEN                                               
            DSQ=AIMAG(Z)**2+(REAL(Z)-REAL(ZPOS))**2                             
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            PHASE(IF)=PHASE(IF)+CON*ATAN2(-AIMAG(Z),REAL(ZPOS)-REAL(Z))         
C                                                                               
C         Real pole.                                                            
          ELSEIF (TYPE.EQ.PZTRP) THEN                                           
            DSQ=AIMAG(Z)**2+(REAL(Z)-REAL(ZPOS))**2                             
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            PHASE(IF)=PHASE(IF)-CON*ATAN2(-AIMAG(Z),REAL(ZPOS)-REAL(Z))         
C                                                                               
C         Complex zero.                                                         
          ELSEIF (TYPE.EQ.PZTCZ) THEN                                           
            DSQ=(REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)-AIMAG(ZPOS))**2               
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(AIMAG(ZPOS)-AIMAG(Z),REAL(ZPOS)-REAL(Z))             
            PHASE(IF)=PHASE(IF)+TEMP                                            
C                                                                               
C         Complex pole.                                                         
          ELSEIF (TYPE.EQ.PZTCP) THEN                                           
            DSQ=(REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)-AIMAG(ZPOS))**2               
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(AIMAG(ZPOS)-AIMAG(Z),REAL(ZPOS)-REAL(Z))             
            PHASE(IF)=PHASE(IF)-TEMP                                            
C                                                                               
C         Complex conjugate zero pair.                                          
          ELSEIF (TYPE.EQ.PZT2Z) THEN                                           
            DSQ=(REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)-AIMAG(ZPOS))**2               
            DSQ=DSQ*((REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)+AIMAG(ZPOS))**2)         
            AMPSQ(IF)=AMPSQ(IF)+10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(AIMAG(ZPOS)-AIMAG(Z),REAL(ZPOS)-REAL(Z))             
            PHASE(IF)=PHASE(IF)+TEMP                                            
            TEMP=CON*ATAN2(-(AIMAG(ZPOS)+AIMAG(Z)),REAL(ZPOS)-REAL(Z))          
            PHASE(IF)=PHASE(IF)+TEMP                                            
C                                                                               
C         Complex conjugate pole pair.                                          
          ELSEIF (TYPE.EQ.PZT2P) THEN                                           
            DSQ=(REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)-AIMAG(ZPOS))**2               
            DSQ=DSQ*((REAL(Z)-REAL(ZPOS))**2+(AIMAG(Z)+AIMAG(ZPOS))**2)         
            AMPSQ(IF)=AMPSQ(IF)-10.0*ALOG10(DSQ)                                
            TEMP=CON*ATAN2(AIMAG(ZPOS)-AIMAG(Z),REAL(ZPOS)-REAL(Z))             
            PHASE(IF)=PHASE(IF)-TEMP                                            
            TEMP=CON*ATAN2(-(AIMAG(ZPOS)+AIMAG(Z)),REAL(ZPOS)-REAL(Z))          
            PHASE(IF)=PHASE(IF)-TEMP                                            
C                                                                               
C         Type unknown.                                                         
          ELSE                                                                  
            GOTO 9004                                                           
C                                                                               
          ENDIF                                                                 
C                                                                               
C       Add on magnitude and phase introduced by gain of singularity.           
          LGAIN=20.0*ALOG10(CABS(ZGAIN))                                        
          PHGAIN=CON*ATAN2(AIMAG(ZGAIN),REAL(ZGAIN))                            
        PHASE(IF)=PHASE(IF)+PHGAIN                                              
        AMPSQ(IF)=AMPSQ(IF)+LGAIN                                               
C                                                                               
5000    CONTINUE                                                                
5100  CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
C                                                                               
C     (4) CALL OF ERROR ROUTINES                                                
C     **************************                                               
C                                                                               
C     Type of FSPEC not s or z.                                                 
9001  CONTINUE                                                                  
      IER=6                                                                     
      GOTO 9999                                                                 
C                                                                               
C     NPOINTS less than one.                                                    
9002  CONTINUE                                                                  
      IER=7                                                                     
      GOTO 9999                                                                 
C                                                                               
C     F1 < F2                                                                   
9003  CONTINUE                                                                  
      IER=8                                                                     
      GOTO 9999                                                                 
C                                                                               
C     Singularity type unknown.                                                 
9004  CONTINUE                                                                  
      IER=9                                                                     
      GOTO 9999                                                                 
C                                                                               
C     FSPEC set up incorrectly - FSPEC(PSN) or FSPEC(PZN) LE 0                  
9005  CONTINUE                                                                  
      IER=10                                                                    
      GOTO 9999                                                                 
C                                                                               
9999  CONTINUE                                                                  
C                                                                               
      CALL UERERR('FILTER',RTN,IER,0,'E')                                       
      CALL UERREP('FILTER',RTN,.TRUE.)                                          
C                                                                               
      RETURN                                                                    
      END                                                                       
