SUBROUTINE FPPJ13
*-- Author :    I. O. Skillicorn      16/11/92
      SUBROUTINE FPPJ13
*     call fpphit                                                                                             
**: FPPJ13 40000 IS. New linking code.                                                                        
**----------------------------------------------------------------------                                      
**: FPPJ13 30907 RP. Farm changes.                                                                            
**----------------------------------------------------------------------                                      
C     JOIN 2 PLANAR MODULES  - 1+3                                                                            
       COMMON/FGMIOS/                                                   
*    PLANAR GEOMETRY                                                                                          
     + ZPP(140),C(140),S(140),ITYPE(140),WZERO(140),WSPACE,             
*                                                                                                             
*    RADIAL GEOMETRY                                                                                          
     + ZP(36),PHW(36),WS(36)                                            
                                                                        
*KEEP,FPJPAR.                                                                                                 
      COMMON/FPJPAR/                                                    
     +      RRCUT1, RRCUT2, RRCUT3,                                     
     +      PLCC3, PLCC12, PLCC23, PLCC13                               
*KEND.                                                                                                        
                                                                        
      COMMON/FTPS3/NS(3),SPAR(4,50,3),IPT(12,50,3),IPLA(12,50,3),       
     1  SGN(12,50,3),YYS(12,50,3),YYF(12,50,3)                          
      COMMON/FTPLNK/LINK,CHIL(100),LNK(3,100),IUS(100,3)                
      COMMON/FTPPBK/NPP,IPP(36,100),CHPP(100),LP(3,100)                 
      COMMON/FPPFIT/PSSS(100),PISS(100),RSSS(100),RISS(100)             
                                                                        
      COMMON/FKLOC/KLOC(100)                                            
      COMMON/FEVSAT/IEVSAT                                              
      COMMON/fsegtp/iseg(100,3)                                         
                                                                        
      DIMENSION PSS(100),PIS(100),RSS(100),RIS(100)                     
      DIMENSION PPP(50),XX(36),YY(36),ZZ(36),WP(36),WPP(36)             
      DIMENSION CH(100)                                                 
      PI2=6.283185307                                                   
*     RRCUT=10.0                                                                                              
      LINK=0                                                            
      LINKO=LINK                                                        
      DO 100 I=1,NS(1)                                                  
      IF(IUS(I,1).NE.0)GOTO100                                          
      Z1=ZPP(6)                                                         
*     X1=SPAR(3,I,1)*Z1+SPAR(4,I,1)                                                                           
*     Y1=SPAR(1,I,1)*Z1+SPAR(2,I,1)                                                                           
      DO 300 K=1,NS(3)                                                  
      IF(IUS(K,3).NE.0)GOTO300                                          
      Z3=ZPP(30)                                                        
*     X3=SPAR(3,K,3)*Z3+SPAR(4,K,3)                                                                           
*     Y3=SPAR(1,K,3)*Z3+SPAR(2,K,3)                                                                           
      ZM=0.5*(Z1+Z3)                                                    
      X1M=SPAR(3,I,1)*ZM+SPAR(4,I,1)                                    
      Y1M=SPAR(1,I,1)*ZM+SPAR(2,I,1)                                    
      X3M=SPAR(3,K,3)*ZM+SPAR(4,K,3)                                    
      Y3M=SPAR(1,K,3)*ZM+SPAR(2,K,3)                                    
      RR=    ((X1M-X3M)**2+(Y1M-Y3M)**2)                                
      CALL SHS(903,0,RR)                                                                               
      RRCUT=RRCUT1                                                      
      If(   iseg(i,1).gt.1                                              
     + .or. iseg(k,3).gt.1  )RRCUT=RRCUT2                               
      If(   iseg(i,1).gt.2                                              
     + .or. iseg(k,3).gt.2  )RRCUT=RRCUT3                               
      IF(RR.GT.RRCUT)GOTO300                                            
      call fpphit(1,3,i,k,iflag)                                        
      if(iflag.eq.1)goto300                                             
      IC=0                                                              
      DO 400 L=1,3                                                      
      IF(L.EQ.2)GOTO400                                                 
      IF(L.EQ.1)II=I                                                    
C      IF(L.EQ.2)II=J                                                                                         
      IF(L.EQ.3)II=K                                                    
      DO 410 LL=1,12                                                    
      IP=(L-1)*12+LL                                                    
      ZA=ZPP(IP)                                                        
      XF=SPAR(3,II,L)*ZA+SPAR(4,II,L)                                   
      YF=SPAR(1,II,L)*ZA+SPAR(2,II,L)                                   
      IC=IC+1                                                           
      XX(IC)=ZA                                                         
      YY(IC)=ATAN2(YF,XF)                                               
      ZZ(IC)=SQRT(XF**2+YF**2)                                          
      IF(YY(IC).LT.0.0)YY(IC)=YY(IC)+PI2                                
      WP(IC)=1./(0.10/ZZ(IC))                                           
      WPP(IC)=1.0                                                       
 410  CONTINUE                                                          
 400  CONTINUE                                                          
      IF(IC .GT. 1)THEN                                                 
      DO600 JJ=2,IC                                                     
      DP = YY(JJ)-YY(JJ-1)                                              
      IF(DP.GT.0.) THEN                                                 
      IF(ABS(DP).GT.ABS(DP-PI2))YY(JJ)=YY(JJ)-PI2                       
      ELSE                                                              
      IF(ABS(DP).GT.ABS(DP+PI2))YY(JJ)=YY(JJ)+PI2                       
      ENDIF                                                             
 600  CONTINUE                                                          
      ENDIF                                                             
C     FIT PHI-Z                                                                                               
      CALL FTLFTW(XX,YY,WP,IC,0,2,PS,PI,D1,D2,D3,D4)
C     FIT R-Z                                                                                                 
      CALL FTLFTW(XX,ZZ,WPP,IC,0,2,RS,RI,D1,D2,D3,D4)
      CHIP=0.                                                           
      DO 700 JJ=1,IC                                                    
      CHIP=CHIP+(YY(JJ)-PS*XX(JJ)-PI)**2*WP(JJ)**2                      
*     DIFF=YY(JJ)-PS*XX(JJ)-PI                                                                                
*     PRINT1002,JJ,(YY(JJ)-PS*XX(JJ)-PI)**2*WP(JJ)**2,YY(JJ),DIFF                                             
 1002 FORMAT(' ',I3,F10.2,2F10.4)                                       
 700  CONTINUE                                                          
      CHIP=CHIP/24.                                                     
      CALL FPCHI(1,3,I,K,CHID)
      CALL SHS(530,0,CHIP)                                                                             
      CALL SHS(574,0,CHID)                                                                             
*     Remove Links with poor Chisq...  18/11/93                                                               
      IF(CHID.GT.100.)GOTO300                                           
*                                                                                                             
      LINK=LINK+1                                                       
      IF(LINK.GT.100)LINK=100                                           
      PSS(LINK)=PS                                                      
      PIS(LINK)=PI                                                      
      RSS(LINK)=RS                                                      
      RIS(LINK)=RI                                                      
      CHIL(LINK)=CHIP                                                   
      CH(LINK)=CHID                                                     
      LNK(1,LINK)=I                                                     
      LNK(2,LINK)=0                                                     
      LNK(3,LINK)=K                                                     
*      WRITE(*,*)' I,K,LINK,CHIP,CHID ',I,K,LINK,CHIP,CHID                                                    
C     MOMENTUM FROM ANGLE CHANGE                                                                              
*     T1=ATAN2(SPAR(1,I,1),SPAR(3,I,1))                                                                       
*     IF(T1.LT.0.)T1=T1+PI2                                                                                   
*     T2=ATAN2(SPAR(1,K,3),SPAR(3,K,3))                                                                       
*     IF(T2.LT.0.)T2=T2+PI2                                                                                   
*     TANT=(SQRT(X3**2+Y3**2)-SQRT(X1**2+Y1**2))/(Z3-Z1)                                                      
*     SL=SQRT((X1-X3)**2+(Y1-Y3)**2)                                                                          
*     THET=ATAN(TANT)                                                                                         
*     DIFF=T2-T1                                                                                              
*     IF(DIFF.LT.-3.14)DIFF=DIFF+PI2                                                                          
*     IF(DIFF.GT.3.14)DIFF=DIFF-PI2                                                                           
*     PP=1000.                                                                                                
*     IF(DIFF.NE.0.0)PP=-0.0002998*12.*SL/(DIFF*SIN(THET))                                                    
*     PPP(LINK)=PP                                                                                            
C     CALL SHS(590,0,1./PP)                                                                                   
C     CALL SHS(591,0,1./PP)                                                                                   
C     CALL SHS(592,0,PP)                                                                                      
*     II=LINK                                                                                                 
*     PRINT 1001,II,LNK(1,II),LNK(2,II),LNK(3,II),CHIL(II)                                                    
 1001 FORMAT(' PJ13I',I3,2X,3I3,F10.2)                                  
 300  CONTINUE                                                          
 200  CONTINUE                                                          
 100  CONTINUE                                                          
                                                                        
C     NEW COMPARE SECTION                                                                                     
      LLL=LINK                                                          
      DO 500 LOOP=1,LLL                                                 
      CHB=100000.                                                       
      KB=0                                                              
C     SELECT BEST                                                                                             
      DO 510 K=1,LLL                                                    
      IF(CH(K).LT.0.0)GOTO 510                                          
      IF(LNK(1,K).EQ.0)GOTO 510                                         
*     IF(LNK(2,K).EQ.0)GOTO 510                                                                               
      IF(LNK(3,K).EQ.0)GOTO 510                                         
*      WRITE(*,*)' K ,CHI  ',K,CH(K),LNK(1,K),LNK(2,K),LNK(3,K)                                               
      IF(CH(K).LT.CHB)THEN                                              
      CHB=CH(K)                                                         
      KB=K                                                              
      ENDIF                                                             
 510  CONTINUE                                                          
*      WRITE(*,*)' KB,CHIB ',KB,CHB                                                                           
      IF(KB.EQ.0)GOTO 545                                               
C     COMPARE BEST WITH REMAINDER                                                                             
      DO 520 K=1,LLL                                                    
      IF(K.EQ.KB)GOTO520                                                
      IF(LNK(1,K).EQ.0)GOTO520                                          
*     IF(LNK(2,K).EQ.0)GOTO520                                                                                
      IF(LNK(3,K).EQ.0)GOTO520                                          
      IF(LNK(1,K).EQ.LNK(1,KB))GOTO530                                  
*     IF(LNK(2,K).EQ.LNK(2,KB))GOTO530                                                                        
      IF(LNK(3,K).EQ.LNK(3,KB))GOTO530                                  
      GOTO 520                                                          
C     REMOVE LINK                                                                                             
 530  LNK(1,K)=0                                                        
      LNK(2,K)=0                                                        
      LNK(3,K)=0                                                        
*      WRITE(*,*)' REMOVE ',K                                                                                 
 520  CONTINUE                                                          
C     COMPARE FINISHED , MARK BEST SEGMENT USED                                                               
      CH(KB)=-CH(KB)                                                    
 500  CONTINUE                                                          
C     RESET CHI WHEN COMPARE FINISHED                                                                         
 545  DO 540 LOOP=1,LLL                                                 
      IF(CH(LOOP).LT.0.0)CH(LOOP)=-CH(LOOP)                             
 540  CONTINUE                                                          
C                                                                                                             
C                                                                                                             
      DO 550 I=LINKO+1,LINK                                             
      IF(LNK(1,I).EQ.0)GOTO550                                          
C     SET USED FLAG                                                                                           
      IUS(LNK(1,I),1)=1                                                 
C     IUS(LNK(2,I),2)=1                                                                                       
      IUS(LNK(3,I),3)=1                                                 
      L1=LNK(1,I)                                                       
C     L2=LNK(2,I)                                                                                             
      L3=LNK(3,I)                                                       
C     PRINT 1000,I,(IPT(III,L1,1),III=1,12)                                                                   
C    1,(IPT(III,L3,3),III=1,12)                                                                               
C    1,CHIL(I),PPP(I)                                                                                         
 1000 FORMAT(' PJN13',I3,1X,2(1X,12I2),F6.2,F7.2)                       
      NPP=NPP+1                                                         
      IF(NPP.GT.100) THEN                                               
        NPP = 100                                                       
        IEVSAT = 1                                                      
      ENDIF                                                             
      KLOC(NPP) = 4                                                     
                                                                        
      LP(1,NPP)=L1                                                      
      LP(2,NPP)=00                                                      
      LP(3,NPP)=L3                                                      
      PSSS(NPP)=PSS(I)                                                  
      PISS(NPP)=PIS(I)                                                  
      RSSS(NPP)=RSS(I)                                                  
      RISS(NPP)=RIS(I)                                                  
      DO 551 KK=1,12                                                    
      IPP(KK,NPP)=IPT(KK,L1,1)                                          
C     IPP(KK+12,NPP)=IPT(KK,L2,2)                                                                             
      IPP(KK+24,NPP)=IPT(KK,L3,3)                                       
 551  CONTINUE                                                          
      CHPP(NPP)=CH(I)                                                   
      PROD=PROB(CH(I)*24.,24)                                           
      CALL SHS(578,0,PROD )                                                                            
      CALL SHS(550,0,8.001)                                                                            
      CALL SHS(550,0,10.001)                                                                           
 550  CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
*