C     PROGRAM SAS TO CALCULATE HAUPTMAN's SAS TRIPLE's ESTIMATE                 
C     CALCULATE THE A VALUE AND OMEGA FOR EACH 3-PHASE INVARIANT                
C     Author: David Langs langs@algol.hwi.buffalo.edu 
C
C     NF = NUMBER OF BIGGEST E-VALUES TO BE PHASED                              
      PARAMETER (NF=4000)                                                       
      INTEGER*2 IH,IK,IL,JK0,JL0,IP1,IP2,ITEM(4)                                
      DOUBLE PRECISION FMH, DELTAH, FMK, DELTAK, FML, DELTAL                    
      INTEGER SER(3),FORM(3)                                                    
      COMMON /ATOMS/ NTYPE, ANUM(6)                                             
      COMMON /AVAL/ NA, AVAL(8), OMEGA(8)                                       
      COMMON /PSCAT/ FMH(6), DELTAH(6), FMK(6), DELTAK(6),                      
     *               FML(6), DELTAL(6)                                          
      DIMENSION F(6,NF),DLT(6,NF),E(2,NF),PH(NF)                                
      DIMENSION CEL(6),DFC(9,6),DT(2,6)                                         
      PI1 = ACOS(-1.0)                                                          
      PI2 = PI1+PI1                                                             
      write(*,11)                                                               
   11 format(' Minimum A-value for triples estimate : ',$)                      
      read(*,*)THR,SCALE                                                        
      if(SCALE.eq.0.0) scale=1.                                                 
      t0 = secnds(0.0)                                                          
      OPEN(1,FILE='sas.dat',STATUS='OLD',FORM='FORMATTED')                      
C     READ CELL CONSTANTS, ANGLES IN DEGREES                                    
      READ(1,*)CEL                                                              
      RAD=PI1/180.                                                              
      CA=COS(CEL(4)*RAD)                                                        
      CB=COS(CEL(5)*RAD)                                                        
      CG=COS(CEL(6)*RAD)                                                        
      SA=SIN(CEL(4)*RAD)                                                        
      SB=SIN(CEL(5)*RAD)                                                        
      SG=SIN(CEL(6)*RAD)                                                        
      V=CEL(1)*CEL(2)*CEL(3)*SQRT(1.-CA**2 -CB**2 -CG**2 +2*CA*CB*CG)           
      AS=CEL(2)*CEL(3)*SA/V                                                     
      BS=CEL(1)*CEL(3)*SB/V                                                     
      CS=CEL(1)*CEL(2)*SG/V                                                     
      A1=2.0*as*bs*(CA*CB-CG)/(SA*SB)                                           
      A2=2.0*as*cs*(CA*CG-CB)/(SA*SG)                                           
      A3=2.0*bs*cs*(CB*CG-CA)/(SB*SG)                                           
C     READ NUMBER OF CHEMICAL ATOM TYPES, AND NUMBER OF ATOMS OF EACH           
      READ(1,*)NTYPE                                                            
      READ(1,*) (ANUM(I),I=1,NTYPE)                                             
C     READ FORM FACTOR COEFFICIENTS FOR EACH ATOM TYPE                          
      DO 4 I=1,NTYPE                                                            
    4 READ(1,*) (DFC(J,I),J=1,9)                                                
      READ(1,*) (DT(1,I),I=1,NTYPE)                                             
      READ(1,*) (DT(2,I),I=1,NTYPE)                                             
      CLOSE(1)                                                                  
      open(1,file='SAS.hkl',status='old',form='formatted')                      
      open(2,file='trip.sas',status='unknown',form='formatted')                 
      DO 10 N=1,4000                                                            
      READ(1,*,END=321) IH,IK,IL,(E(I,N),I=1,2),PH(N)                           
      stl = (ih*as)**2 +(ik*bs)**2 +(il*cs)**2                                  
      s = 0.25*(stl +a1*ih*ik +a2*ih*il + a3*ik*il)                             
      DO 5 L=1,NTYPE                                                            
      F0=DFC(1,L)*EXP(-DFC(2,L)*S)+DFC(3,L)*EXP(-DFC(4,L)*S)+                   
     *DFC(5,L)*EXP(-DFC(6,L)*S)+DFC(7,L)*EXP(-DFC(8,L)*S)+DFC(9,L)              
      SS=DT(2,L)                                                                
      CC=F0+DT(1,L)                                                             
      DLT(L,N)=ATAN2(SS,CC)                                                     
      F(L,N)=SQRT(SS*SS+CC*CC)                                                  
    5 CONTINUE                                                                  
      if(N.lt.5) write(*,2)ih,ik,il,(f(i,n),i=1,ntype),                         
     *(dlt(i,N),i=1,ntype),s                                                    
   2  FORMAT(6x,3i3,6f9.4/6f9.6,3f5.2)                                          
      PH(N)=AMOD(PH(N)+PI2,PI2)                                                 
      IF(PH(N).GT.PI1) PH(N)=PH(N)-PI2                                          
  10  CONTINUE                                                                  
  321 NFX=N-1                                                                   
      FORM(1)=1                                                                 
      N1=0                                                                      
      open(3,file='trip.tmp',form='unformatted')                                
      DO 3 itrip=1,9999999                                                      
      read(3,end=322)item                                                       
      it=item(4)                                                                
      IH=item(1)                                                                
      IK=abs(item(2))                                                           
      IL=abs(item(3))                                                           
      FORM(2)=item(2)/ik                                                        
      FORM(3)=item(3)/il                                                        
      DO 40 J=1,NTYPE                                                           
      FMH(J)=F(J,IH)                                                            
      FMK(J)=F(J,IK)                                                            
      FML(J)=F(J,IL)                                                            
      DELTAH(J)=DLT(J,IH)                                                       
      DELTAK(J)=DLT(J,IK)                                                       
  40  DELTAL(J)=DLT(J,IL)                                                       
      PSI=PH(IH)                                                                
      JIT=IT.AND.15                                                             
      PSI = PI2+PH(IH)+FORM(2)*PH(IK)+FORM(3)*PH(IL) + JIT*PI1/6.               
      PSI = MOD(PSI,PI2)                                                        
      IF(PSI.gt.PI1) PSI=PSI-PI2                                                
      SER(1)=IH                                                                 
      SER(2)=IK                                                                 
      SER(3)=IL                                                                 
      CALL AVALUE(SER,FORM,E,NFX)                                               
      JK0=JK                                                                    
      JL0=JL                                                                    
      AV = SCALE*(AVAL(1)+AVAL(2))/2.                                           
      IF(AV.LT.THR) GOTO 3                                                      
      N1=N1+1                                                                   
      IP1=1000.0*PSI                                                            
      IP2=500.0*(OMEGA(1)-OMEGA(2))                                             
      WRITE(2,1)AV,IH,ITEM(2),ITEM(3),IT,IP2,IP1                                
   1  format(f7.3,i5,2i6,3i6)                                                   
   3  continue                                                                  
 322  T1 = secnds(t0)                                                           
      WRITE(*,100) NFX,N1,T1                                                    
 100  FORMAT(I7,' REFLECTIONS',I9,' ACCEPTED INVARIANTS. T= ',F5.0)             
      STOP                                                                      
      END                                                                       
C     ************************************************************              
      SUBROUTINE AVALUE (SERIAL,FORM,E,NF)                                      
C                                                                               
      PARAMETER ( MXLC = 32000, MXRF = 5000 )                                   
C                                                                               
      INTEGER*2 MRP, MS, MP, MPSER, ISINTL, NSER                                
      INTEGER SERIAL(3), FORM(3)                                                
C                                                                               
      DOUBLE PRECISION FMH, DELTAH, FMK, DELTAK, FML, DELTAL                    
      DOUBLE PRECISION FMH2, DELH2, FMK2, DELK2, FML2, DELL2, FMHKL             
      DOUBLE PRECISION A2H, A2K, A2L, CH, SH, CK, SK, CL, SL                    
      DOUBLE PRECISION X, X1, X2, X3, XI, XI1, XI2, XI3                         
      DOUBLE PRECISION CHKL, SHKL, CHBKL, SHBKL, CHKBL, SHKBL,                  
     *                 CHKLB, SHKLB, AAA, ANG                                   
      DOUBLE PRECISION GAMMA, GA0, GA1, GA2, GA3                                
      DOUBLE PRECISION SIGMA, SI0, SI1, SI2, SI3                                
      DOUBLE PRECISION Z, Z0, Z1, Z2, Z3, ZETA, ZETA0, ZETA1, ZETA2,            
     *                 ZETA3                                                    
      DOUBLE PRECISION C0, S0, CPART, SPART, B, XXX                             
      DOUBLE PRECISION R, R1, R2, R3, RB, R1B, R2B, R3B, R123, R1B23,           
     *                 R12B3, R123B, R1B2B3, R1B23B, R12B3B, R1B2B3B            
      DOUBLE PRECISION TAU, T1, T2, T3, T12, T13, T23, T123                     
C                                                                               
      DIMENSION FMH2(15), DELH2(15), FMK2(15), DELK2(15), FML2(15),             
     *          DELL2(15), FMHKL(15)                                            
      DIMENSION X(3), XI(3), TAU(3)                                             
      DIMENSION Z(0:3), ZETA(0:3), GAMMA(0:3), SIGMA(0:3)                       
      DIMENSION R(3), RB(3),E(2,NF)                                             
C                                                                               
      COMMON /ATOMS/ NTYPE, ANUM(6)                                             
      COMMON /AVAL/ NA, AVAL(8), OMEGA(8)                                       
      COMMON /PSCAT/ FMH(6), DELTAH(6), FMK(6), DELTAK(6),                      
     *               FML(6), DELTAL(6)                                          
      COMMON NSORT,MAXTR,MTRP,ISTR1,ISTR2,IOPT,NGROUP,I2D                       
      COMMON NBV,NBVT,NLIST,NSPAC,EMIN,AMIN,EMINR,AMINR,AMABS                   
      COMMON N3R,TITLE(27),ET2,ET3,ETT,AMAXTR                                   
      COMMON MHKL(MXRF),MRP(MXRF),MS(MXRF),MP(MXRF),MPSER(MXRF)                 
      COMMON ISINTL(MXRF),NSER(MXLC)                                            
C                                                                               
      EQUIVALENCE (X(1),X1), (X(2),X2), (X(3),X3)                               
      EQUIVALENCE (XI(1),XI1), (XI(2),XI2), (XI(3),XI3)                         
      EQUIVALENCE (GAMMA(0),GA0), (GAMMA(1),GA1), (GAMMA(2),GA2),               
     *            (GAMMA(3),GA3)                                                
      EQUIVALENCE (SIGMA(0),SI0), (SIGMA(1),SI1), (SIGMA(2),SI2),               
     *            (SIGMA(3),SI3)                                                
      EQUIVALENCE (Z(0),Z0), (Z(1),Z1), (Z(2),Z2), (Z(3),Z3)                    
      EQUIVALENCE (ZETA(0),ZETA0), (ZETA(1),ZETA1),                             
     *            (ZETA(2),ZETA2), (ZETA(3),ZETA3)                              
      EQUIVALENCE (TAU(1),T1), (TAU(2),T2), (TAU(3),T3)                         
      EQUIVALENCE (R(1),R1), (R(2),R2), (R(3),R3)                               
      EQUIVALENCE (RB(1),R1B), (RB(2),R2B), (RB(3),R3B)                         
C                                                                               
      DATA PI2/6.2831853/                                                       
C                                                                               
C Set up scattering factor expresions which remain constant.                    
C                                                                               
      DO 200 I=1,NTYPE                                                          
         FMH2(I) = ANUM(I)*FMH(I)*FMH(I)                                        
         FMK2(I) = ANUM(I)*FMK(I)*FMK(I)                                        
         FML2(I) = ANUM(I)*FML(I)*FML(I)                                        
         DELH2(I) = 2.*DELTAH(I)                                                
         DELK2(I) = 2.*DELTAK(I)                                                
         DELL2(I) = 2.*DELTAL(I)                                                
         FMHKL(I) = ANUM(I)*FMH(I)*FMK(I)*FML(I)                                
  200 CONTINUE                                                                  
C                                                                               
C Calculate "ALPHA sub 2H", "ALPHA sub K", etc.                                 
C                                                                               
      A2H = 0.                                                                  
      A2K = 0.                                                                  
      A2L = 0.                                                                  
      DO 205 I=1,NTYPE                                                          
         A2H = A2H + FMH2(I)                                                    
         A2K = A2K + FMK2(I)                                                    
         A2L = A2L + FML2(I)                                                    
  205 CONTINUE                                                                  
C                                                                               
C Calculate "C sub H", "S sub K", etc..                                         
C                                                                               
      CH = 0.                                                                   
      SH = 0.                                                                   
      CK = 0.                                                                   
      SK = 0.                                                                   
      CL = 0.                                                                   
      SL = 0.                                                                   
      DO 210 I=1,NTYPE                                                          
         CH = CH + FMH2(I)*COS(DELH2(I))                                        
         SH = SH + FMH2(I)*SIN(DELH2(I))                                        
         CK = CK + FMK2(I)*COS(DELK2(I))                                        
         SK = SK + FMK2(I)*SIN(DELK2(I))                                        
         CL = CL + FML2(I)*COS(DELL2(I))                                        
         SL = SL + FML2(I)*SIN(DELL2(I))                                        
  210 CONTINUE                                                                  
      CH = CH / A2H                                                             
      SH = SH / A2H                                                             
      CK = CK / A2K                                                             
      SK = SK / A2K                                                             
      CL = CL / A2L                                                             
      SL = SL / A2L                                                             
C                                                                               
C Obatin "X sub 1", "XI sub 1", "X sub 2", "Xi sub 2", etc..                    
C                                                                               
      X1 = SQRT(CH*CH + SH*SH)                                                  
      XI1 = ATAN2(-SH/X1,CH/X1)                                                 
      X2 = SQRT(CK*CK + SK*SK)                                                  
      XI2 = ATAN2(-SK/X2,CK/X2)                                                 
      X3 = SQRT(CL*CL + SL*SL)                                                  
      XI3 = ATAN2(-SL/X3,CL/X3)                                                 
C                                                                               
C Calculate "C sub HKL", S sub HKL", etc..                                      
C                                                                               
      CHKL  = 0.                                                                
      SHKL  = 0.                                                                
      CHBKL = 0.                                                                
      SHBKL = 0.                                                                
      CHKBL = 0.                                                                
      SHKBL = 0.                                                                
      CHKLB = 0.                                                                
      SHKLB = 0.                                                                
      DO 220 I=1,NTYPE                                                          
         ANG   =  DELTAH(I) + DELTAK(I) + DELTAL(I)                             
         CHKL  = CHKL + FMHKL(I)*COS(ANG)                                       
         SHKL  = SHKL + FMHKL(I)*SIN(ANG)                                       
         ANG   = -DELTAH(I) + DELTAK(I) + DELTAL(I)                             
         CHBKL = CHBKL + FMHKL(I)*COS(ANG)                                      
         SHBKL = SHBKL + FMHKL(I)*SIN(ANG)                                      
         ANG   =  DELTAH(I) - DELTAK(I) + DELTAL(I)                             
         CHKBL = CHKBL + FMHKL(I)*COS(ANG)                                      
         SHKBL = SHKBL + FMHKL(I)*SIN(ANG)                                      
         ANG   =  DELTAH(I) + DELTAK(I) - DELTAL(I)                             
         CHKLB = CHKLB + FMHKL(I)*COS(ANG)                                      
         SHKLB = SHKLB + FMHKL(I)*SIN(ANG)                                      
  220 CONTINUE                                                                  
      AAA = 1. / SQRT(A2H*A2K*A2L)                                              
      CHKL  = CHKL *AAA                                                         
      SHKL  = SHKL *AAA                                                         
      CHBKL = CHBKL*AAA                                                         
      SHBKL = SHBKL*AAA                                                         
      CHKBL = CHKBL*AAA                                                         
      SHKBL = SHKBL*AAA                                                         
      CHKLB = CHKLB*AAA                                                         
      SHKLB = SHKLB*AAA                                                         
C                                                                               
C Calculate the Gamma's and Sigma's.                                            
C                                                                               
      GA0= CHKL *(1. - (CH*CK*CL - CH*SK*SL - SH*CK*SL - SH*SK*CL) ) +          
     1     SHKL *      (SH*SK*SL - SH*CK*CL - CH*SK*CL - CH*CK*SL)   +          
     2     CHBKL*(-CH + (CK*CL - SK*SL) ) +                                     
     3     SHBKL*( SH + (CK*SL + SK*CL) ) +                                     
     4     CHKBL*(-CK + (CH*CL - SH*SL) ) +                                     
     5     SHKBL*( SK + (CH*SL + SH*CL) ) +                                     
     6     CHKLB*(-CL + (CH*CK - SH*SK) ) +                                     
     7     SHKLB*( SL + (CH*SK + SH*CK) )                                       
C                                                                               
      SI0= CHKL *      (SH*SK*SL - SH*CK*CL - CH*SK*CL - CH*CK*SL)   +          
     1     SHKL *(1. + (CH*CK*CL - CH*SK*SL - SH*CK*SL - SH*SK*CL) ) +          
     2     CHBKL*(-SH + (CK*SL + SK*CL) ) +                                     
     3     SHBKL*(-CH - (CK*CL - SK*SL) ) +                                     
     4     CHKBL*(-SK + (CH*SL + SH*CL) ) +                                     
     5     SHKBL*(-CK - (CH*CL - SH*SL) ) +                                     
     6     CHKLB*(-SL + (CH*SK + SH*CK) ) +                                     
     7     SHKLB*(-CL - (CH*CK - SH*SK) )                                       
C                                                                               
      GA1= CHKL *(-CH + (CK*CL - SK*SL) ) +                                     
     1     SHKL *(-SH + (CK*SL + SK*CL) ) +                                     
     2     CHBKL*(1. - (CH*CK*CL - CH*SK*SL + SH*CK*SL + SH*SK*CL) ) +          
     3     SHBKL*(    -(SH*SK*SL - SH*CK*CL + CH*SK*CL + CH*CK*SL) ) +          
     4     CHKBL*(-CL + (CH*CK + SH*SK) ) +                                     
     5     SHKBL*(-SL - (CH*SK - SH*CK) ) +                                     
     6     CHKLB*(-CK + (CH*CL + SH*SL) ) +                                     
     7     SHKLB*(-SK - (CH*SL - SH*CL) )                                       
C                                                                               
      SI1= CHKL *( SH + (CK*SL + SK*CL) ) +                                     
     1     SHKL *(-CH - (CK*CL - SK*SL) ) +                                     
     2     CHBKL*(    -(SH*SK*SL - SH*CK*CL + CH*SK*CL + CH*CK*SL) ) +          
     3     SHBKL*(1. + (CH*CK*CL - CH*SK*SL + SH*CK*SL + SH*SK*CL) ) +          
     4     CHKBL*(-SL + (CH*SK - SH*CK) ) +                                     
     5     SHKBL*( CL + (CH*CK + SH*SK) ) +                                     
     6     CHKLB*(-SK + (CH*SL - SH*CL) ) +                                     
     7     SHKLB*( CK + (CH*CL + SH*SL) )                                       
C                                                                               
      GA2= CHKL *(-CK + (CH*CL - SH*SL) ) +                                     
     1     SHKL *(-SK + (CH*SL + SH*CL) ) +                                     
     2     CHBKL*(-CL + (CK*CH + SK*SH) ) +                                     
     3     SHBKL*(-SL - (CK*SH - SK*CH) ) +                                     
     4     CHKBL*(1. - (CH*CK*CL + CH*SK*SL - SH*CK*SL + SH*SK*CL) ) +          
     5     SHKBL*(    -(SH*SK*SL + SH*CK*CL - CH*SK*CL + CH*CK*SL) ) +          
     6     CHKLB*(-CH + (CK*CL + SK*SL) ) +                                     
     7     SHKLB*(-SH - (CK*SL - SK*CL) )                                       
C                                                                               
      SI2= CHKL *( SK + (CH*SL + SH*CL) ) +                                     
     1     SHKL *(-CK - (CH*CL - SH*SL) ) +                                     
     2     CHBKL*(-SL + (CK*SH - SK*CH) ) +                                     
     3     SHBKL*( CL + (CK*CH + SK*SH) ) +                                     
     4     CHKBL*(    -(SH*SK*SL + SH*CK*CL - CH*SK*CL + CH*CK*SL) ) +          
     5     SHKBL*(1. + (CH*CK*CL + CH*SK*SL - SH*CK*SL + SH*SK*CL) ) +          
     6     CHKLB*(-SH + (CK*SL - SK*CL) ) +                                     
     7     SHKLB*( CH + (CK*CL + SK*SL) )                                       
C                                                                               
      GA3= CHKL *(-CL + (CH*CK - SH*SK) ) +                                     
     1     SHKL *(-SL + (CH*SK + SH*CK) ) +                                     
     2     CHBKL*(-CK + (CL*CH + SL*SH) ) +                                     
     3     SHBKL*(-SK - (CL*SH - SL*CH) ) +                                     
     4     CHKBL*(-CH + (CL*CK + SL*SK) ) +                                     
     5     SHKBL*(-SH - (CL*SK - SL*CK) ) +                                     
     6     CHKLB*(1. - (CH*CK*CL + CH*SK*SL + SH*CK*SL - SH*SK*CL) ) +          
     7     SHKLB*(    -(SH*SK*SL + SH*CK*CL + CH*SK*CL - CH*CK*SL) )            
C                                                                               
      SI3= CHKL *( SL + (CH*SK + SH*CK) ) +                                     
     1     SHKL *(-CL - (CH*CK - SH*SK) ) +                                     
     2     CHBKL*(-SK + (CL*SH - SL*CH) ) +                                     
     3     SHBKL*( CK + (CL*CH + SL*SH) ) +                                     
     4     CHKBL*(-SH + (CL*SK - SL*CK) ) +                                     
     5     SHKBL*( CH + (CL*CK + SL*SK) ) +                                     
     6     CHKLB*(    -(SH*SK*SL + SH*CK*CL + CH*SK*CL - CH*CK*SL) ) +          
     7     SHKLB*(1. + (CH*CK*CL + CH*SK*SL + SH*CK*SL - SH*SK*CL) )            
C                                                                               
C Obtain the "Z sub j" and "Sigma sub j".                                       
C                                                                               
      DO 230 J=0,3                                                              
         Z(J) = SQRT(GAMMA(J)*GAMMA(J) + SIGMA(J)*SIGMA(J))                     
         ZETA(J) = ATAN2(SIGMA(J)/Z(J),GAMMA(J)/Z(J))                           
  230 CONTINUE                                                                  
C                                                                               
C Obtain "Tau sub j".                                                           
C                                                                               
      DO 235 J=1,3                                                              
         IF (FORM(J).GT.0) THEN                                                 
            R (J) = E(1,SERIAL(J))                                              
            RB(J) = E(2,SERIAL(J))                                              
         ELSE                                                                   
            R (J) = E(2,SERIAL(J))                                              
            RB(J) = E(1,SERIAL(J))                                              
         END IF                                                                 
         ARG = 2.*R(J)*RB(J)*X(J) / (1. - X(J)*X(J))                            
      TAJ = 0.0                                                                 
      AA = ARG*ARG                                                              
      IF(ARG.GT.4.0) TAJ = 1.0 - 0.5/ARG - 0.16/AA                              
      IF(ARG.LT.1.0) TAJ = 0.50287*ARG -.015708*AA -0.041078*ARG*AA             
      IF(TAJ.EQ.0.0) TAJ = 0.57745*ARG -0.13869*AA +0.012091*ARG*AA             
         TAU(J) = TAJ                                                           
  235 CONTINUE                                                                  
C                                                                               
C Set up constant parts of equations.                                           
C                                                                               
      T12  = T1*T2                                                              
      T13  = T1*T3                                                              
      T23  = T2*T3                                                              
      T123 = T1*T23                                                             
C                                                                               
      R123    = R1*R2*R3                                                        
      R1B23   = R1B*R2*R3                                                       
      R12B3   = R1*R2B*R3                                                       
      R123B   = R1*R2*R3B                                                       
      R1B2B3  = R1B*R2B*R3                                                      
      R1B23B  = R1B*R2*R3B                                                      
      R12B3B  = R1*R2B*R3B                                                      
      R1B2B3B = R1B*R2B*R3B                                                     
C                                                                               
C XXX is used in calculating A values.                                          
C                                                                               
      XXX = 1. / ( (1. - X1*X1) * (1. - X2*X2) * (1. - X3*X3) )                 
C                                                                               
C Calculate A values and OMEGA's for the eight invariants.                      
C                                                                               
      C0 = Z0*( R123           * COS(ZETA0) +                                   
     1          R1B2B3B * T123 * COS(XI1 + XI2 + XI3 + ZETA0) ) +               
     2     Z1*( R1B23   * T1   * COS(XI1 - ZETA1) +                             
     3          R12B3B  * T23  * COS(XI2 + XI3 + ZETA1) ) +                     
     4     Z2*( R12B3   * T2   * COS(XI2 - ZETA2) +                             
     5          R1B23B  * T13  * COS(XI1 + XI3 + ZETA2) ) +                     
     6     Z3*( R123B   * T3   * COS(XI3 - ZETA3) +                             
     7          R1B2B3  * T12  * COS(XI1 + XI2 + ZETA3) )                       
C                                                                               
      S0 = Z0*( R123           * SIN(ZETA0) -                                   
     1          R1B2B3B * T123 * SIN(XI1 + XI2 + XI3 + ZETA0) ) -               
     2     Z1*( R1B23   * T1   * SIN(XI1 - ZETA1) +                             
     3          R12B3B  * T23  * SIN(XI2 + XI3 + ZETA1) ) -                     
     4     Z2*( R12B3   * T2   * SIN(XI2 - ZETA2) +                             
     5          R1B23B  * T13  * SIN(XI1 + XI3 + ZETA2) ) -                     
     6     Z3*( R123B   * T3   * SIN(XI3 - ZETA3) +                             
     7          R1B2B3  * T12  * SIN(XI1 + XI2 + ZETA3) )                       
         CPART = C0                                                             
         SPART = S0                                                             
         B = SQRT(CPART*CPART + SPART*SPART)                                    
         AVAL (1) = 2.*B*XXX                                                    
         OMEGA(1) = ATAN2(SPART/B,CPART/B)                                      
C                                                                               
      C0 = Z0*( R1B2B3B        * COS(ZETA0) +                                   
     1          R123    * T123 * COS(XI1 + XI2 + XI3 + ZETA0) ) +               
     2     Z1*( R12B3B  * T1   * COS(XI1 - ZETA1) +                             
     3          R1B23   * T23  * COS(XI2 + XI3 + ZETA1) ) +                     
     4     Z2*( R1B23B  * T2   * COS(XI2 - ZETA2) +                             
     5          R12B3   * T13  * COS(XI1 + XI3 + ZETA2) ) +                     
     6     Z3*( R1B2B3  * T3   * COS(XI3 - ZETA3) +                             
     7          R123B   * T12  * COS(XI1 + XI2 + ZETA3) )                       
C                                                                               
      S0 = Z0*( R1B2B3B        * SIN(ZETA0) -                                   
     1          R123    * T123 * SIN(XI1 + XI2 + XI3 + ZETA0) ) -               
     2     Z1*( R12B3B  * T1   * SIN(XI1 - ZETA1) +                             
     3          R1B23   * T23  * SIN(XI2 + XI3 + ZETA1) ) -                     
     4     Z2*( R1B23B  * T2   * SIN(XI2 - ZETA2) +                             
     5          R12B3   * T13  * SIN(XI1 + XI3 + ZETA2) ) -                     
     6     Z3*( R1B2B3  * T3   * SIN(XI3 - ZETA3) +                             
     7          R123B   * T12  * SIN(XI1 + XI2 + ZETA3) )                       
         CPART = C0                                                             
         SPART = S0                                                             
         B = SQRT(CPART*CPART + SPART*SPART)                                    
         AVAL (2) = 2.*B*XXX                                                    
         OMEGA(2) = ATAN2(SPART/B,CPART/B)                                      
      RETURN                                                                    
      END                                                                       
                                                                                
