      SUBROUTINE BPF(SSPECI,SSPECO,FL,FU)                                       
C                                                                               
C  MODULE         BPF                                                           
C  VERSION        1.0                                                           
C  ROUTINE NO     70                                                            
C  LAST MODIFIED  21/12/81                                                      
C  AUTHOR         DAVE POWELL                                                   
C                                                                               
C  THIS SUBROUTINE TAKES A PROTOTYPE LOW PASS FILTER WITH A                     
C  BREAKPOINT AT ANGULAR FREQUENCY 1RAD/SEC AND CONVERTS TO                     
C  A BANDPASS FILTER WITH LOW FREQUENCY BREAK POINT AT f=FL                     
C  AND UPPER BREAKPOINT f=FU.                                                   
C                                                                               
C  PARAMETERS                                                                   
C  **********                                                                   
C                                                                               
C  SSPECI   R*4A  INPUT S-PLANE SPECIFICATION                                   
C  SSPECO   R*4A  OUTPUT S-PLANE SPECIFICATION                                  
C  FL       R*4   LOWER BREAK FREQUENCY                                         
C  FU       R*4   UPPER BREAK FREQUENCY                                         
C                                                                               
C  THE LENGTH OF SSPECI AND SSPECO ARE INDETERMINANT                            
      REAL*4 FL,FU,SSPECI(9999),SSPECO(9999)                                    
C                                                                               
C  GLOBAL VARIABLES                                                             
C  ****************                                                             
C                                                                               
$INCLUDE SIGPRO10.CMN (NLIST)                                                   
$INCLUDE SPLANE10.CMN (NLIST)                                                   
C                                                                               
C  LOCAL VARIABLES                                                              
C  ***************                                                              
C                                                                               
      INTEGER*4 RTN,NSI,INFNO,IPSI,ITS,INOS,TYPE,NSO                            
      REAL*4 PI,WL,WU,DIFF,PROD,PROD4,RPOS                                      
      COMPLEX CPROD4,CTEM,CTEMSQ,SPOS,SGAIN                                     
      PARAMETER (RTN=70)                                                        
C                                                                               
C  (1) CHECK PARAMETERS                                                         
C  ********************                                                         
C                                                                               
      IF (SSPECI(PSX).NE.PSV) GOTO 9001                                         
      IF (FL.LT.0.0.OR.FU.LT.0.0) GOTO 9007                                     
      IF (FL.GE.FU) GOTO 9002                                                   
C                                                                               
C  (2) INITIALISE                                                               
C  **************                                                               
C                                                                               
      NSI=INT(SSPECI(PSN))                                                      
      IF (NSI.LT.0) GOTO 9003                                                   
      PI=4.0*ATAN(1.0)                                                          
      WL=2.0*PI*FL                                                              
      WU=2.0*PI*FU                                                              
      DIFF=WU-WL                                                                
      PROD=WU*WL                                                                
      PROD4=4.0*PROD                                                            
      INFNO=0                                                                   
      INOS=0                                                                    
C                                                                               
C  (3) DETERMINE LENGTH OF SSPECO                                               
C  ******************************                                               
C                                                                               
      IF (NSI.LE.0) GOTO 3100                                                   
      DO 3000 IS=1,NSI                                                          
          IPSI=PSH+1+PSE*(IS-1)                                                 
C         Check sequence number of singularity.                                 
          IF (INT(SSPECI(IPSI+PSM)).NE.IS) GOTO 9004                            
C         Increment counter according to S-Plane singularity type.              
          ITS=INT(SSPECI(IPSI+PST))                                             
C                                                                               
C         Real zero.                                                            
          IF (ITS.EQ.PSTRZ) THEN                                                
              INFNO=INFNO-1                                                     
              RPOS=SSPECI(IPSI+PSR)                                             
              IF ((RPOS*DIFF)**2.GT.PROD4) THEN                                 
                  INOS=INOS+2                                                   
              ELSE                                                              
                  INOS=INOS+1                                                   
              ENDIF                                                             
C                                                                               
C         Real Pole.                                                            
          ELSEIF (ITS.EQ.PSTRP) THEN                                            
              INFNO=INFNO+1                                                     
              RPOS=SSPECI(IPSI+PSR)                                             
              IF ((RPOS*DIFF)**2.GT.PROD4) THEN                                 
                  INOS=INOS+2                                                   
              ELSE                                                              
                  INOS=INOS+1                                                   
              ENDIF                                                             
C                                                                               
C         Complex Zero.                                                         
          ELSEIF (ITS.EQ.PSTCZ) THEN                                            
              INFNO=INFNO-1                                                     
              INOS=INOS+2                                                       
C                                                                               
C         Complex Pole.                                                         
          ELSEIF (ITS.EQ.PSTCP) THEN                                            
              INFNO=INFNO+1                                                     
              INOS=INOS+2                                                       
C                                                                               
C         Complex Conjugate Zero Pair.                                          
          ELSEIF (ITS.EQ.PST2Z) THEN                                            
              INFNO=INFNO-2                                                     
              INOS=INOS+2                                                       
C                                                                               
C         Complex Conjugate Pole Pair.                                          
          ELSEIF (ITS.EQ.PST2P) THEN                                            
              INFNO=INFNO+2                                                     
              INOS=INOS+2                                                       
C                                                                               
C         Singularity type not known.                                           
          ELSE                                                                  
              GOTO 9006                                                         
          ENDIF                                                                 
3000  CONTINUE                                                                  
3100  CONTINUE                                                                  
C                                                                               
C  (4) INITIALISE HEADER OF SSPECO                                              
C  *******************************                                              
C                                                                               
C     Find number of singularities in output (SSPECO)                           
      IF ((INFNO/2)*2.EQ.INFNO) THEN                                            
          NSO=IABS(INFNO)/2+INOS                                                
      ELSE                                                                      
          NSO=(IABS(INFNO)+1)/2+INOS                                            
      ENDIF                                                                     
      SSPECO(PSN)=NSO                                                           
C                                                                               
C     Set type of output array (SSPECO).                                        
      SSPECO(PSX)=PSV                                                           
C                                                                               
C     Check output (SSPECO) is long enough.                                     
      IF (PSH+PSE*NSO.GT.INT(SSPECO(PSL))) THEN                                 
          GOTO 9005                                                             
      ENDIF                                                                     
C                                                                               
C     Change max used length if necessary.                                      
      IF (PSH+PSE*NSO.GT.INT(SSPECO(PSU))) THEN                                 
          SSPECO(PSU)=PSH+PSE*NSO                                               
      ENDIF                                                                     
C                                                                               
C  (5) CALCULATION FOR EACH SINGULARITY-OUTPUT INTO SSPECO                      
C  *******************************************************                      
C                                                                               
      IF (NSO.LE.0) THEN                                                        
          GOTO 5100                                                             
      ENDIF                                                                     
      ISO=0                                                                     
C                                                                               
C     Sweep over singularities.                                                 
      DO 5000 ISI=1,NSI                                                         
          IPSI=PSH+1+PSE*(ISI-1)                                                
          SPOS=CMPLX(SSPECI(IPSI+PSR),SSPECI(IPSI+PSI))                         
          TYPE=INT(SSPECI(IPSI+PST))                                            
          SGAIN=CMPLX(SSPECI(IPSI+PSG),SSPECI(IPSI+PSJ))                        
          TEM=REAL(SPOS)*DIFF                                                   
          TEMSQ=TEM**2                                                          
          CTEM=SPOS*DIFF                                                        
          CTEMSQ=CTEM**2                                                        
          CPROD4=CMPLX(PROD4,0.0)                                               
C                                                                               
C         Real Zero.                                                            
          IF (TYPE.EQ.PSTRZ) THEN                                               
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              IF (TEMSQ.GT.PROD4) THEN                                          
                  SSPECO(IPSO+PSM)=ISO                                          
                  SSPECO(IPSO+PST)=PSTRZ                                        
                  SSPECO(IPSO+PSG)=REAL(SGAIN)/DIFF                             
                  SSPECO(IPSO+PSJ)=AIMAG(SGAIN)/DIFF                            
                  SSPECO(IPSO+PSR)=(TEM+SQRT(TEMSQ-PROD4))/2.0                  
                  SSPECO(IPSO+PSI)=0.0                                          
                  ISO=ISO+1                                                     
                  IPSO=PSH+1+PSE*(ISO-1)                                        
                  SSPECO(IPSO+PSM)=ISO                                          
                  SSPECO(IPSO+PST)=PSTRZ                                        
                  SSPECO(IPSO+PSG)=1.0                                          
                  SSPECO(IPSO+PSJ)=0.0                                          
                  SSPECO(IPSO+PSR)=(TEM-SQRT(TEMSQ-PROD4))/2.0                  
                  SSPECO(IPSO+PSI)=0.0                                          
              ELSE                                                              
                  SSPECO(IPSO+PSM)=ISO                                          
                  SSPECO(IPSO+PST)=PST2Z                                        
                  SSPECO(IPSO+PSG)=SSPECI(IPSI+PSG)/DIFF                        
                  SSPECO(IPSO+PSJ)=0.0                                          
                  SSPECO(IPSO+PSR)=TEM/2.0                                      
                  SSPECO(IPSO+PSI)=(SQRT(PROD4-TEMSQ))/2.0                      
              ENDIF                                                             
C                                                                               
C         Real Pole.                                                            
          ELSEIF (TYPE.EQ.PSTRP) THEN                                           
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              IF (TEMSQ.GT.PROD4) THEN                                          
                  SSPECO(IPSO+PSM)=ISO                                          
                  SSPECO(IPSO+PST)=PSTRP                                        
                  SSPECO(IPSO+PSG)=REAL(SGAIN)*DIFF                             
                  SSPECO(IPSO+PSJ)=AIMAG(SGAIN)*DIFF                            
                  SSPECO(IPSO+PSR)=(TEM+SQRT(TEMSQ-PROD4))/2.0                  
                  SSPECO(IPSO+PSI)=0.0                                          
                  ISO=ISO+1                                                     
                  IPSO=PSH+1+PSE*(ISO-1)                                        
                  SSPECO(IPSO+PST)=PSTRP                                        
                  SSPECO(IPSO+PSG)=1.0                                          
                  SSPECO(IPSO+PSJ)=0.0                                          
                  SSPECO(IPSO+PSR)=(TEM-SQRT(TEMSQ-PROD4))/2.0                  
                  SSPECO(IPSO+PSI)=0.0                                          
              ELSE                                                              
                  SSPECO(IPSO+PSM)=ISO                                          
                  SSPECO(IPSO+PST)=PST2P                                        
                  SSPECO(IPSO+PSG)=REAL(SGAIN)*DIFF                             
                  SSPECO(IPSO+PSJ)=AIMAG(SGAIN)*DIFF                            
                  SSPECO(IPSO+PSR)=TEM/2.0                                      
                  SSPECO(IPSO+PSI)=(SQRT(PROD4-TEMSQ))/2.0                      
              ENDIF                                                             
C                                                                               
C         Complex Zero.                                                         
          ELSEIF (TYPE.EQ.PSTCZ) THEN                                           
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PSTCZ                                            
              SSPECO(IPSO+PSG)=REAL(SGAIN)/DIFF                                 
              SSPECO(IPSO+PSJ)=AIMAG(SGAIN)/DIFF                                
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS+CSQRT(TEMSQ-CPROD4))/2.0          
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS+CSQRT(TEMSQ-CPROD4))/2.0         
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PSTCZ                                            
              SSPECO(IPSO+PSG)=1.0                                              
              SSPECO(IPSO+PSJ)=0.0                                              
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0        
C                                                                               
C         Complex Pole.                                                         
          ELSEIF (TYPE.EQ.PSTCP) THEN                                           
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PSTCP                                            
              SSPECO(IPSO+PSG)=REAL(SGAIN)*DIFF                                 
              SSPECO(IPSO+PSJ)=AIMAG(SGAIN)*DIFF                                
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS+CSQRT(TEMSQ-CPROD4))/2.0          
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS+CSQRT(TEMSQ-CPROD4))/2.0         
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PSTCP                                            
              SSPECO(IPSO+PSG)=1.0                                              
              SSPECO(IPSO+PSJ)=0.0                                              
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0        
C                                                                               
C         Complex Conjugate Zero Pair.                                          
          ELSEIF (TYPE.EQ.PST2Z) THEN                                           
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PST2Z                                            
              SSPECO(IPSO+PSG)=REAL(SGAIN)/(DIFF**2)                            
              SSPECO(IPSO+PSJ)=AIMAG(SGAIN)/(DIFF**2)                           
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS+CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS+CSQRT(CTEMSQ-CPROD4))/2.0        
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PSG)=1.0                                              
              SSPECO(IPSO+PSJ)=0.0                                              
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0        
C                                                                               
C         Complex Conjugate Pole Pair.                                          
          ELSEIF (TYPE.EQ.PST2P) THEN                                           
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PST2P                                            
              SSPECO(IPSO+PSG)=REAL(SGAIN)*(DIFF**2)                            
              SSPECO(IPSO+PSJ)=AIMAG(SGAIN)*(DIFF**2)                           
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS+CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS+CSQRT(CTEMSQ-CPROD4))/2.0        
              ISO=ISO+1                                                         
              IPSO=PSH+1+PSE*(ISO-1)                                            
              SSPECO(IPSO+PSM)=ISO                                              
              SSPECO(IPSO+PST)=PST2P                                            
              SSPECO(IPSO+PSG)=1.0                                              
              SSPECO(IPSO+PSJ)=0.0                                              
              SSPECO(IPSO+PSR)=REAL(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0         
              SSPECO(IPSO+PSI)=AIMAG(DIFF*SPOS-CSQRT(CTEMSQ-CPROD4))/2.0        
C                                                                               
          ENDIF                                                                 
5000  CONTINUE                                                                  
5100  CONTINUE                                                                  
C                                                                               
C  (6) APPEND PAIRS OF POLES/ZEROES AT S=0                                      
C  ***************************************                                      
C                                                                               
      IF (IABS(INFNO)/2.LE.0) GOTO 6100                                         
      DO 6000 IS=1,IABS(INFNO)/2                                                
          ISO=ISO+1                                                             
          IPSO=PSH+1+PSE*(ISO-1)                                                
          IF (INFNO.GT.0) THEN                                                  
              SSPECO(IPSO+PST)=PST2Z                                            
          ELSE                                                                  
              SSPECO(IPSO+PST)=PST2P                                            
          ENDIF                                                                 
          SSPECO(IPSO+PSM)=ISO                                                  
          SSPECO(IPSO+PSG)=1.0                                                  
          SSPECO(IPSO+PSJ)=0.0                                                  
          SSPECO(IPSO+PSR)=0.0                                                  
          SSPECO(IPSO+PSI)=0.0                                                  
6000  CONTINUE                                                                  
6100  CONTINUE                                                                  
C                                                                               
C  (7) APEND SINGLE REMAINING POLE/ZERO IF IT EXISTS                            
C  *************************************************                            
C                                                                               
      IF ((INFNO/2)*2.EQ.INFNO) GOTO 7100                                       
      ISO=ISO+1                                                                 
      IPSO=PSH+1+PSE*(ISO-1)                                                    
      IF (INFNO.GT.0) THEN                                                      
          SSPECO(IPSO+PST)=PSTRZ                                                
      ELSE                                                                      
          SSPECO(IPSO+PST)=PSTRP                                                
      ENDIF                                                                     
      SSPECO(IPSO+PSM)=ISO                                                      
      SSPECO(IPSO+PSG)=1.0                                                      
      SSPECO(IPSO+PSJ)=0.0                                                      
      SSPECO(IPSO+PSR)=0.0                                                      
      SSPECO(IPSO+PSI)=0.0                                                      
7100  CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
C                                                                               
C  (8) ERRORS                                                                   
C  **********                                                                   
C                                                                               
C     Input array format not s-plane.                                           
9001  CONTINUE                                                                  
      IER=20                                                                    
      GOTO 9999                                                                 
C                                                                               
C     Lower breakpoint greater than upper breakpoint.                           
9002  CONTINUE                                                                  
      IER=21                                                                    
      GOTO 9999                                                                 
C                                                                               
C     SSPECI(PSN) is negative.                                                  
9003  CONTINUE                                                                  
      IER=22                                                                    
      GOTO 9999                                                                 
C                                                                               
C     Error in sequence number in input array.                                  
9004  CONTINUE                                                                  
      IER=23                                                                    
      GOTO 9999                                                                 
C                                                                               
C     Output array is not long enough.                                          
9005  CONTINUE                                                                  
      IER=24                                                                    
      GOTO 9999                                                                 
C                                                                               
C     Unknown type of singularity.                                              
9006  CONTINUE                                                                  
      IER=25                                                                    
      GOTO 9999                                                                 
C                                                                               
C     FL or FU is less than zero.                                               
9007  CONTINUE                                                                  
      IER=26                                                                    
      GOTO 9999                                                                 
C                                                                               
9999  CONTINUE                                                                  
C                                                                               
      CALL UERERR('FILTER',RTN,IER,0,'E')                                       
      CALL UERREP('FILTER',RTN,.TRUE.)                                          
C                                                                               
      RETURN                                                                    
      END                                                                       
