C  PROGRAM PATMAP   PROGRAMMED BY BOB SNYDER  5/15/67  FOR I.B.M. 1620          
C  PATMAP IS A GENERAL PROGRAM FOR CALCULATING AND PLOTTING A TWO DIMENSIONAL   
C  PATTERSON MAP PROJECTED ALONG ONE OF THE THREE PRINCIPAL CRYSTALLOGRAPHIC    
C  AXIS.  THE TWO PRINCIPAL LIMITATIONS OF PATMAP ARE                           
C    1.THE EQUATION PROGRAAMED IS FOR THE GENERAL TRICLINIC CASE.               
C    THEREFORE ALL REFLECTIONS MUST BE USED NOT JUST INDEPENDANT REFLECTIONS    
C    IF ONLY INDEPENDANT REFLECTIONS  WERE TAKEN THEN THE VARIABLE SUBROUTINE   
C    FUNCT SHOULD BE CHANGED TO CONTAIN THE EQUATION SPECIFIC TO YOUR           
C    SPACE GROUP( VOL1 OF INTERNATIONAL TABLES)                                 
C    2.  THE TWO DIMENSIONAL PROJECTION MAY ONLY BE ALONG ONE OF THE THREE      
C    PRINCIPAL CRYSTALLOGRAPHIC AXIS                                            
C  NOTE ON THE 1620 THIS PROGRAM IS VERY SLOW TAKING 3 MIN TO PRINT ONE LINE    
C  OF OUTPUT WHEN 150 REFLECTIONS ARE USED.  THEREFORE SPECIFIED SHOULD         
C  NOT BE MUCH SMALLER THAN .01                                                 
C*******************************************************************************
C                    DATA SET UP                                                
C*******************************************************************************
C*******************************************************************************
C  CARD1  DIRECT CELL CONSTANTS  (6F10.6)                                       
C        COL 1-10  A IN ANGSTROMS                                               
C        COL11-20  B *    *                                                     
C        COL21-30  C *    *                                                     
C        COL31-40  COS ALPHA                                                    
C        COL41-50  COS BETA                                                     
C        COL51-60  COS GAMMA                                                    
C*******************************************************************************
C  CARD 2   LIMITS   (6F10.2)                                                   
C        COL 1-10  UPPER LIMIT OF X                                             
C        COL11-20  LOWER LIMIT OF X                                             
C        COL21-30  UPPER LIMIT OF Y                                             
C        COL31-40  LOWER LIMIT OF Y                                             
C        COL41-50  UPPER LIMIT OF Z                                             
C        COL51-60  LOWER LIMIT OF Z                                             
C  NOTE  THE LIMITS FOR THE PROJECTION AXIS MUST BE ZERO                        
C*******************************************************************************
C  CARD 3     (3F5.0,I5)                                                        
C        COL 1- 5  THE INTERVAL DESIRED BETWEEN EACH X VALUE                    
C        COL 6-10  THE INTERVAL DESIRED BETWEEN EACH Y VALUE                    
C        COL11-15 THE INTERVAL DESIRED BETWEEN EACH Z VALUE                     
C        COL16-20  NEG  =-1  ALL VALUES WILL BE PRINTED ON OUTPUT               
C                       =0  ALL NEGATIVE VALUES WILL BE OMITTED                 
C                       =+NUMBER  ONLY VALUES GREATER THAN NUMBER WILL BE PRINTE
C*******************************************************************************
C  CARD 4 (AND ON)   DATA CARDS  (3I9,F9.2)                                     
C        COL 1- 9  H                                                            
C        COL10-18   K                                                           
C        COL19-27   L                                                           
C        COL28-36   F SQUARED                                                   
C*******************************************************************************
C  FOLLOW THE LAST CARD WITH A CARD PUNCHED ONLY IN COL 1                       
C  LAST CARD (I1)                                                               
C        COL1   =1 MEANS READ 4 NEW SET UP CARDS AND A NEW BATCH OF DATA        
C               =2  MEANS CALL EXIT                                             
C*******************************************************************************
C*******************************************************************************
      COMMON N(250),M(250),FSQ(250),PUV(21),IPUV(21),AN(250),AM(250)            
      IN=105                                                                    
      ITP=108                                                                   
  998 READ(IN,3)A,B,C,CA,CB,CG                                                  
    3 FORMAT(6F10.6)                                                            
      WRITE(ITP,4)A,B,C,CA,CB,CG                                                
    4 FORMAT(15H0CELL CONSTANTS//3H A=F10.6,3H B=F10.6,3H C=F10.6,              
     111H COS ALPHA=F10.6, 10H COS BETA=F10.6, 11H COS GAMMA=F10.6)             
      READ(IN,40)UX,BX,UY,BY,UZ,BZ                                              
   40 FORMAT(6F10.2)                                                            
      READ(IN,61)TNTX,TNTY,TNTZ,NEG                                             
   61 FORMAT(3F5.0,I5)                                                          
      IF(NEG)90,90,91                                                           
   90 WRITE(ITP,86)                                                             
   86 FORMAT(46H0NEGATIVE VALUES WILL NOT BE PRINTED ON OUTPUT)                 
      GO TO 92                                                                  
   91 WRITE(ITP,87)NEG                                                          
   87 FORMAT(17H0ALL VALUES BELOW  I5, 24H WILL BE PRINTED AS ZERO)             
   92 WRITE(ITP,44)UX,BX,TNTX                                                   
   44 FORMAT(7H0LIMITS//12H X GOES FROMF5.2,3H TOF5.2,12H IN STEPS OFF5.        
     13)                                                                        
      WRITE(ITP,45)UY,BY,TNTY                                                   
   45 FORMAT(12H Y GOES FROMF5.2,3H TOF5.2,12H IN STEPS OFF5.3)                 
      WRITE(ITP,46)UZ,BZ,TNTZ                                                   
   46 FORMAT(12H Z GOES FROMF5.2,3H TOF5.2,12H IN STEPS OFF5.3)                 
      I=1                                                                       
      IF(UX)6,5,6                                                               
    5 IF(BX)6,7,6                                                               
    7 UX=UY                                                                     
      BX=BY                                                                     
      UY=UZ                                                                     
      BY=BZ                                                                     
      TNTX=TNTY                                                                 
      TNTY=TNTZ                                                                 
      WRITE(ITP,8)                                                              
    8 FORMAT(30H0X ON OUTPUT MEANS Y DIRECTION)                                 
      WRITE(ITP,9)                                                              
    9 FORMAT(30H0Y ON OUTPUT MEANS Z DIRECTION)                                 
   10 READ(IN,11)NUM,N(I),M(I),FSQ(I)                                           
       IF(FSQ(I).LT.0.0)GO TO 10                                                
   11 FORMAT(I1,8XI9,I9,F9.2)                                                   
      I=I+1                                                                     
      IF(NUM)10,10,20                                                           
    6 IF(UY)13,12,13                                                            
   12 IF(BY)13,14,13                                                            
   14 UY=UZ                                                                     
      BY=BZ                                                                     
      TNTY=TNTZ                                                                 
      WRITE(ITP,15)                                                             
   15 FORMAT(30H0Y ON OUTPUT MEANS Z DIRECTION)                                 
   16 READ(IN,17)NUM,N(I),M(I),FSQ(I)                                           
      IF(FSQ(I).LT.0.0)GO TO 16                                                 
   17 FORMAT(I1,I8,9XI9,F9.2)                                                   
      I=I+1                                                                     
      IF(NUM)16,16,20                                                           
   13 READ(IN,18)NUM,N(I),M(I),FSQ(I)                                           
      IF(FSQ(I).LT.0.0)GO TO 13                                                 
   18 FORMAT(I1,I8,I9,9XF9.2)                                                   
      I=I+1                                                                     
      IF(NUM)13,13,20                                                           
   20 I=I-2                                                                     
      WRITE(ITP,49)I                                                            
   49 FORMAT(52H0THE NUMBER OF REFLECTIONS TO BE USED IN SUMMATION =I3)         
      VOL= A*B*C* SQRT ( 1.0+2.0 *CA*CB*CG-CA**2 -CB**2 -CG**2)                 
      WRITE(ITP,60)VOL                                                          
   60 FORMAT(30H THE VOLUME OF THE UNIT CELL =F10.2)                            
      WRITE(ITP,21)                                                             
   21 FORMAT(///)                                                               
      ITAG=0                                                                    
      PI=3.1415927                                                              
      U=BX                                                                      
      V=UY                                                                      
      IU=0                                                                      
      RAD=2*PI                                                                  
      POO=0.0                                                                   
      DO 25 J=1,I                                                               
      AN(J)=N(J)                                                                
      AM(J)=M(J)                                                                
   25 POO=POO+FSQ(J)                                                            
      ANORM=999.*VOL/POO                                                        
   30 DO 42 J=1,21                                                              
      IPUV(J)=0.0                                                               
   42 PUV(J)=0.0                                                                
      CALL FUNCT(IU,I,RAD,U,V,UX,UY,BX,BY,TNTX,TNTY)                            
      UN=U                                                                      
      IU=IU-1                                                                   
   62 DO28 IN=1,IU                                                              
      PUV(IN)= PUV(IN)*ANORM/VOL                                                
   28 IPUV(IN)=PUV(IN)                                                          
      CALL OUTPT (V,UX,BY,BX,UY,IU,TNTX,TNTY,NEG,ITAG)                          
      IF(V-BY)31,31,80                                                          
   80 U=BX                                                                      
      V=V-TNTY                                                                  
      GOTO 30                                                                   
   31 IF(UN-UX)71,70,70                                                         
   71 BX=U                                                                      
      V=UY                                                                      
      WRITE(ITP,21)                                                             
      GO TO 30                                                                  
   70 IF(NUM-1)998,998,999                                                      
  999 CALL EXIT                                                                 
      END                                                                       
      SUBROUTINE FUNCT(IU,I,RAD,U,V,UX,UY,BX,BY,TNTX,TNTY)                      
C  FUNCT IS A VARIABLE SUBROUTINE WHICH SHOULD BE CHANGED TO INCLUDE THE        
C  SPACE GROUP SPECIFIC PATTERSON EQUATION.  THE EQUATION PRESENTLY PROGRAMMED  
C  IS FOR THE GENERAL TRICLINIC CASE P(U,V)=SUM OVERN*SUM OVER M OF             
C  COS(2*PI(N*U+M*V))  WHERE N AND M ARE THE TWO INDICES TO BE VARIED           
C  NOTE   LEAVE THE VOLUME TERM OF THE EQUATION TO BE PROGRAMMED                
      COMMON N(250),M(250),FSQ(250),PUV(21),IPUV(21),AN(250),AM(250)            
      DO26 IU=1,21                                                              
      DO27 IA=1,I                                                               
       PUV(IU)= PUV(IU) + FSQ(IA)* COS (RAD*(AN(IA)*U+AM(IA)*V))                
   27 CONTINUE                                                                  
      IF(U-UX)63,63,62                                                          
   63 U=U+TNTX                                                                  
   26 CONTINUE                                                                  
   62 RETURN                                                                    
      END                                                                       
      SUBROUTINE OUTPT (V,UX,BY,BX,UY,IU,TNTX,TNTY,NEG,ITAG)                    
      COMMON N(250),M(250),FSQ(250),PUV(21),IPUV(21),AN(250),AM(250)            
      ITP=108                                                                   
      ITAG=ITAG+1                                                               
      IF(V-UY)51,50,51                                                          
   50 WRITE(ITP,52)                                                             
   52 FORMAT(61H1                                                 **** X        
     1 ****)                                                                    
      WRITE(ITP,53)                                                             
   53 FORMAT(/)                                                                 
   54 DO55 J=1,IU                                                               
      AJ=J                                                                      
   55 PUV(J)=BX+TNTX*(AJ-1.0)                                                   
      WRITE(ITP,56) (PUV(J),J=1,IU)                                             
   56 FORMAT(10H1       X=21F5.2)                                               
      WRITE(ITP,57)                                                             
   57 FORMAT(60H        ************************************************        
     1****   57H********************************************************        
     2*)                                                                        
      WRITE(ITP,62)                                                             
      WRITE(ITP,62)                                                             
      GO TO 60                                                                  
   51 IF(ITAG-13)60,58,58                                                       
   58 WRITE(ITP,59)                                                             
   59 FORMAT(1H1)                                                               
      ITAG=0                                                                    
      GO TO 54                                                                  
   60 IF(NEG)90,91,91                                                           
   91 DO92 J=1,IU                                                               
      IF(NEG)40,40,41                                                           
   40 IF(IPUV(J))45,45,92                                                       
   45 IPUV(J)=0                                                                 
   41 IF(IPUV(J)-NEG)46,46,92                                                   
   46 IPUV(J)=0                                                                 
   92 CONTINUE                                                                  
   90 WRITE(ITP,61)V,(IPUV(J),J=1,IU)                                           
   61 FORMAT(3H Y=F5.2,2H**21I5,2H**)                                           
      WRITE(ITP,62)                                                             
      WRITE(ITP,62)                                                             
      IF(ITAG)47,70,47                                                          
   47 WRITE(ITP,62)                                                             
   62 FORMAT(10H        **)                                                     
   70 RETURN                                                                    
      END                                                                       
