
                 SUBROUTINE  PDFMAP(IN,IOUT,LINE,NPAGE,IMARK)                 
C                                                                       
C   PROGRAM  FOR  CALCULATIONS  OF PROBABILITY DENSITY  MAPS UP TO SIXTH
C   ORDER CUMULANT PRODUCTS IN CASE OF A GRAM-CHARLIER SERIES EXPANSION.
C   IN CASE OF AN EDGEWORTH SERIES  EXPANSIONS ONLY  COMPUTATIONS  UP TO
C   FOURTH ORDER CUMULANT PRODUCTS ARE AVAILABLE, BUT FOR NORMAL AND FOR
C   SO CALLED EXTENDED MAPS.     A  POINT ARRAY  OUTPUT  CAN BE COMPUTED
C   FOR EACH RANK  SEPERATLY  AND FOR  EACH  COMBINATION OF RANKS OF THE
C   CUMULANTS.                                                          
C   IN ADDITION THE CORRECTLY  NORMALIZED  PROBABILITY  DENSITY FUNCTION
C   OF THE ALPHA-FORMALISM  (BASED ON GRAM-CHARLIER) UP TO  FOURTH ORDER
C   CUMULANT PRODUCTS CAN BE CALCULATED (NOV. 1980)                     
C   PROBABILITY DENSITIES AND POTENTIALS IN ARBITRARY PLANES CAN BE CAL-
C             CULATED SINCE ITSIMPLEMENTATION IN DEC. 1980.             
C                                                                       
C              PROGRAMER    UDO ZUCKER  MAX-PLANCK INSTITUT             
C              FUER  FESTKOERPERFORSCHUNG   7 STUTTGART FRG             
C                          (SEPTEMBER  1979)                            
C              REVISED AND UPDATED             MARCH   1980             
C              REVISED AND ADDED               NOV,DEC 1980             
C              MC-error calculation   W.F.K.   FEB     1990               
C                                                                       
C   INPUT FILES ARE                                                     
C   ---------------                                                     
C                   MAINFILE   OF PROGRAM REFINE            UNIT   ISC1 
C                   HARMO FILE OF PROGRAM REFINE            UNIT IHARMO 
C                   STANDARD    INPUT       UNIT            UNIT     IN 
C                                                                       
C   OUTPUT FILES ARE                                                    
C   ----------------                                                    
C                                                                       
C                   PLOT OUTPUT UNIT FOR ARBITRARY                      
C                   MAPS AND PLOTS                          UNIT S2,S3,S4,
C                                                           S5,S6 
C                   STANDARD   OUTPUT         UNIT          UNIT   IOUT 
C                                                                       
C    DEFINITION OF OTHER INPUT VARIABLES                                
C    -----------------------------------                                
C                                                                       
C    LINE   ARRAY OF DIM. 11.  IT CARRIES INFORMATION ABOUT THAT MAN WHO
C           TAKES CARE ON THE PROGRAMME SYSTEM AT YOUR COMPUTING CENTER.
C    NPAGE  MAXIMUM  NUMBER  OF  LINES  PER PAGE AT THE CURRENT PRINTER.
C                                                                       
      character*25         jdruck,jdr(5)
      DOUBLE PRECISION     JCUM,IAUSG,IFUENF,ISECHS,JNOT                
C      double precision     qp(10000),qm(10000)
      INTEGER              S2,S3,S4,S5,S6,ZRO,ENDD,TITLE,FRMT 
      integer*4            seed
      REAL                 MIROR3,ISKEW                                 
C                                                                       
      DIMENSION   JCUM(5),IAUSG(5),IFUENF(21),ISECHS(28),FRMT(50),IU(50)
     &,LABEL(3), RECELL(6), G(3,3), RMAT(3,3,24),TRANS(3,24),CELLP(3,4),
     &IPARA(1000),FCURVE(9,20),SCALE(10), JSCRF(10),PAR(60),PARNAM(2,60)
     &,IPAR(60), IRPAR(60),EGGY(6),KEXTRF(7), BDF(16),INTER(20),Y(35,4),
     &Z(35,4),EXTRAN(12),ANOM(2,20),EXTIN(6),TAG(2,60),ER(21,6),FR(28,6)
     &,ISOT(60),ISCAT(2,60), OCCA(60),OCCUP(60),SITE(60),XYZ(3,60),BETA(
     &6,60),CR(10,60),DR(15,60), KDEGN(2640),IPDFF(6),XH(3,60),BH(6,60),
     &R1(3), R2(6), HJKLM(15),  X(3),  IP(43),  FX(21), TITLE(18), IPLOT
     &(60),CELL(6), IDAT(5), ISEC(6), ISER(20),IDATA(11),IORIGN(2),ITYPE
     &(12),ALF(60), IOUTP(8),IPRI(34),OUTP(2),JFIVE(21),JSIX(28),LOOP(6)
     &,JTHREE(10),JFOUR(15),LINE(11),E1(3),E2(3),E3(3),MTYP(3),VECTR1(3)
     &,VECTR2(3),POUT(3),XHELP(3)                                       
      dimension cv(40000),bffd(200),jirl(5040),jir(84,60),
     &cl(3570),ranr(85),ranc(85),vcv(84,84),c(84,84),
     &occao(60),xyzo(3,60),betao(6,60),cro(10,60),dro(15,60),ero(21,6),
     &fro(6,28),qs(10000),iminus(10000),iplus(10000),array(10000)
      dimension qp(10000),qm(10000)
C                                                                       
C                                                                       
      COMMON           /PRINT/WINKEL(3),ISELEC,ISECT,IS,IHILF,LSPCH,GAM 
      COMMON           /INVBET/Q(3,3)                                   
      COMMON           /WQ/W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)
      COMMON           /BACH/ARRAY,qp,qm,qs,iplus,iminus                               
      common           /vaco/ cv,bffd
      common           /parms/ isot,occa,xyz,beta,cr,dr,er,fr,jir,ipara
      common           /parmol/ occao,xyzo,betao,cro,dro,ero,fro
      common           /con/ cnstr(1000)
C                                                                       
      EQUIVALENCE      (ARRAY(2641),ISCAT),  
     &(ARRAY(3161),RMAT),(ARRAY(3977),PAR),(ARRAY(4037),PARNAM)  
     &,(ARRAY(4157),IPAR),(ARRAY(4217),IRPAR), (ARRAY(4277),Y), (ARRAY  
     &(4417),Z)                  
      EQUIVALENCE  (WINKEL(1),ALFA),(WINKEL(2),BET),(IOUTP(1),OUTP(1))
      equivalence (jirl(1),jir(1,1))  
C     EQUIVALENCE  (DR(1,41),ER),(DR(7,49),FR)                          
C                                                                       
C----------------- ASSIGNEMENT OF LOGICAL UNIT NUMBERS ---------------- 
C                                                                       
      DATA            ISC1,IHARMO,S2,S3,S4,S5,S6  /10,9,31,32,33,34,35/
C                                                                       
C---------------------------------------------------------------------- 
C                                                                       
      DATA     ISPC,IBNK,TITLE(18),JNOT,IARB/1H0,1H ,4HPDF ,1H ,4HARB / 
      DATA     LEER ,ZRO,NZRO,ENDD              /1H ,3H,4X,3H,I4,1H)/   
      DATA     IDAT /4HNO  ,4HYES ,4HDOUB,4HSING,4HLE  /                
      DATA     KLVIN,MCENTR /1HK,4HCENT/                                
      DATA     ISEC   /4HXYZD, 4HXZYU, 4HYXZU, 4HYZXD, 4HZYXU,4HZXYD/   
      DATA     ISER/4HEDGE,4H-WOR,4HTH-S,4HERIE,4H    ,4HGRAM,4H-CHA,   
     &        4HRLIE,4HR-SE,4HRIE ,4HALPH,4HA-FO,4HRM. ,4HEDGE,4HW.  ,  
     &        4HALPH,4HA-FO,4HRM. ,4HGRAM,4H-CH./                       
      DATA     JCUM/7H  TWO  ,7H THREE ,7H  FOUR ,7H FIVE  ,7H  SIX  /  
      DATA IORIGN,TWOPI3,TWOPI4,TWOPI6/4HZERO,4HATOM,248.050214,1558.   
     &         54546,61528.908/                                         
C                                                                       
C---- IF YOU THINK ABOUT INCONSISTENCY INTERPRET TWOPI2,TWOPI3,.......   
C                                                                       
      DATA     TWOPI2,EINDRI,ZWIDRI/19.7392088,0.33333333,0.66666666/   
      DATA     ITYPE  /4H  FI, 4HRST , 4H SEC, 4HOND , 4H  TH,4HIRD ,   
     &         4H  FO,4HRTH ,4H  FI,4HFTH ,4H  SI,4HXTH /               
      DATA    IPRI/2HX=,2HY=,2HZ=,4HB11=,4HB22=,4HB33=,4HB12=,4HB13=,   
     & 4HB23=,4HC111,4HC222,4HC333,4HC112,4HC122,4HC113,4HC133,4HC223   
     &,4HC233,4HC123,4H1111,4H2222,4H3333,4H1112,4H1113,4H1122,4H1123   
     &,4H1133,4H1222,4H1223,4H1233,4H1333,4H2223,4H2233,4H2333/         
      DATA  IFUENF/7HE11111=,7HE22222=,7HE33333=,7HE11112=,7HE11113=,   
     &             7HE11122=,7HE11123=,7HE11133=,7HE11222=,7HE11223=,   
     &             7HE11233=,7HE11333=,7HE12222=,7HE12223=,7HE12233=,   
     &   7HE12333=,7HE13333=,7HE22223=,7HE22233=,7HE22333=,7HE23333=/   
      DATA  ISECHS/7HF111111,7HF222222,7HF333333,7HF111112,7HF111113,   
     &   7HF111122,7HF111123,7HF111133,7HF111222,7HF111223,7HF111233,   
     &   7HF111333,7HF112222,7HF112223,7HF112233,7HF112333,7HF113333,   
     &   7HF122222,7HF122223,7HF122233,7HF122333,7HF123333,7HF133333,   
     &   7HF222223,7HF222233,7HF222333,7HF223333,7HF233333/             
      DATA LOOP   /3,6,10,15,21,28/                                     
      DATA FRMT     /4H(A1,,4HF4.0,3H,1X,47*1H /                        
      DATA IU       /3H,1X,3H,2X,3H,3X,3H,4X,3H,5X,3H,6X,3H,7X,3H,8X,   
     & 3H,9X,4H,10X,4H,11X,4H,12X,4H,13X,4H,14X,4H,15X,4H,16X,4H,17X,   
     &4H,18X,4H,19X,4H,20X,4H,21X,4H,22X,4H,23X,4H,24X,4H,25X,4H,26X,   
     &4H,27X,4H,28X,4H,29X,4H,30X,4H,31X,4H,32X,4H,33X,4H,34X,4H,35X,   
     &4H,36X,4H,37X,4H,38X,4H,39X,4H,40X,4H,41X,4H,42X,4H,43X,4H,44X,   
     &4H,45X,4H,46X,4H,47X,4H,48X,4H,49X,4H,50X/ 
      data jdr /'                         ',
     &          ' ERROR-PDF (1-sigma) NEG.',
     &          ' ERROR-PDF (1-sigma)     ',
     &          ' 1-sigma SIGNIFICANCE-MAP',
     &          ' N-sigma SIGNIFICANCE-MAP' /                 
C                                                                       
C                                                                       
C   +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
C   +              QUASI DYNAMICALLY CORE ALOCATION                   + 
C   +  MAXIMA MUST BE EQUAL THE DIMENSION OF THE ARRAY NAMED "ARRAY"  + 
C   +                        OF COMMON BACH.                          + 
C   +                                                                 + 
                             MAXIMA=  10000                             
C   +                                                                 + 
C   +       WARNING. MAXIMA MAY NOT LESS THAN 1760 ELEMENTS.          + 
C   +                                                                 + 
C   +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
C                                                                       
C                                                                       
C---- STANDARD OUTPUT ANGLE BETWEEN AXIS IS GAM=90. DEGREES             
C     IEXT=0/1 MEANS NORMAL/EXTENDED   PDF-MAP                          
C                                                                       
      GAM   = 90.                                                       
      IEXT  = 0
      twopi = 2.*acos(-1.)
      seed = 555555
C                                                                       
C---- INDICES FOR STANDARD SELECTION OF A XYZ(DOWN) SECTION             
C                                                                       
      INDEX1= 1                                                         
      INDEX2= 2                                                         
      INDEX3= 3                                                         
C                                                                       
      MIROR3=-1.                                                        
C---- CAUTION  MIROR3 IS REAL VARIABLE. MINUS ONE MEANS Z IS DOWN.      
C                                                                       
C---- READ THE CONTROL CARD.                                            
C                                                                       
      READ(IN,1)(TITLE(L),L=1,17),IPDFF,ISECT,IPLT,LSPC,NPNT,IOR,IX,    
     &     IXI,IXF,IY,IYI,IYF,IZ,IZI,IZF,mc,sifa,isqq,EXXAG,IMINI,IMAXI,
     &     ISELEC,IUPL,MAIN,ICCH,T,KELVIN,PNORM
C---- IF IZ=0 SPECIAL PLANE SHOULB BE CALCULATED.  READ POINTS DE -     
C     FINING PLANE.     
      ISQ = 1                                                
      IF(ISQQ.NE.0)ISQ=0
      IF(IZ.EQ.0) READ(IN,2)MTYP(1),E1,MTYP(2),E2,MTYP(3),E3            
C---- READ ATOM SELECTION SWITCHES.                                     
      READ(IN,3)IPLOT                                                   
    1 FORMAT(17A4/5I1,1X,6I1,10I5,f6.3,i1/,F9.0,3I5,I2,1X,I2,I5,F8.0,
     &A1,F8.0)               
    2 FORMAT(A4,2X,3F8.0/A4,2X,3F8.0/A4,2X,3F8.0)
    3 FORMAT(60I1)                                                      
C                                                                       
C---- WAS THE TEMPERATURE GIVEN BY  T  NOT EQUAL ZERO THE  POTENTIAL    
C     INSTEAD OF THE PDF SHOULD BE CALCULATED.                          
C                                                                       
      IPOTET = 0                                                        
      INORM  = 0                                                        
      IF(ABS(T).LT.1.E-8) GOTO 4                                        
      IPOTET = 1                                                        
      IF(KELVIN.NE.KLVIN) T = T + 273.16                                
      IF(PNORM.LT.1.E-4) INORM = 0                                      
C                                                                       
C---- IF IZ=0 SPECIAL PLANE SHOULD BE CALCULATED.                       
C                                                                       
    4 IPLANE = 0                                                        
      IF(IZ.NE.0) GOTO 128                                              
C---- SET PARAMETERS FOR SPECIAL PLANE TO SENSEFUL VALUES.              
C     INITIALIZE VALUES.                                                
      IPLANE = 1                                                        
      ISECT  = 0                                                        
      ISELEC= 90                                                        
      ITYP1  = 0                                                        
      IOR    = 1                                                        
      IZI    = 0                                                        
      IZF    = 0                                                        
      CELL(3)= 0.                                                       
      XORG   = 0.                                                       
      YORG   = 0.                                                       
      ZORG   = 0.                                                       
C----      JNDEX  = 0   ; NOW INITIALIZED IN MC-LOOP
      IERR   = 0                                                        
C                                                                       
C---- IS ONE OF THE POINTS A CENTER POINT? IF YES, THEN SHUFFLE IT      
C     TO E3.                                                            
      DO 2401 I = 1,3                                                   
      JJ = I                                                            
      IF(MTYP(I).EQ.MCENTR) GOTO 2402                                   
 2401 JJ = 0                                                            
 2402 IF(JJ) 2405,2405,2403                                             
 2403 IF(JJ-3) 2406,2405,2405                                           
 2406 IF(JJ-2) 2409,2407,2407                                           
C---- JJ IS TWO                                                         
 2407 DO 2408 I= 1,3                                                    
      R1(I) = E2(I)                                                     
      E2(I) = E3(I)                                                     
 2408 E3(I) = R1(I)                                                     
      GOTO 2412                                                         
 2409 DO 2410 I = 1,3                                                   
      R1(I) = E1(I)                                                     
      E1(I) = E3(I)                                                     
 2410 E3(I) = R1(I)                                                     
 2412 ITYP1 = 1                                                         
 2405 CONTINUE                                                          
C                                                                       
C---- IF ICCH=0 INCLUDE C(IJK)*C(PQR)*H(IJKPQR) IN CALCULATION OF       
C     AN EDGEWORTH MAP OF FOURTH ORDER.                                 
C     IF ICCH=1  DO NOT INCLUDE THIS EQUATION.                          
C     IF ICCH=2  DO ONLY INCLUDE THIS EQUATION.                         
C                                                                       
  128 ICCH1=0                                                           
      IF(ICCH.EQ.2) ICCH1=1                                             
      IF(ICCH.EQ.2) ICCH =0                                             
C                                                                       
C---- SPECIAL INPUT OR OUTPUT UNIT READ IN ?                                      
C                                                                       
      IF(IUPL.NE.0) S2=IUPL
      write(*,7777) S2
 7777 FORMAT(' S2 = ',I4)
      IF(MAIN)129,130,129                                               
  129 ISC1 = MAIN                                                       
C                                                                       
C---- READ BINARY OUTPUT UNIT 'MAIN' OF THE PROGRAM REFINE.             
C                                                                       
  130 REWIND ISC1                                                       
       READ  (ISC1)                                                     
     &KCYL,LABEL,RECELL,G,VOL,NSYM,NATOM,NSCAT,JCELL,MRANK,NPARA,NSF,DMI
     &N,DMAX,DELMAX,ICENT,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM),((TRANS
     &(I,J),I=1,3),J=1,NSYM),((CELLP(I,J),I=1,3),J=1,JCELL),(IPARA(I),I=
     &1,NPARA),((FCURVE(I,J),I=1,9),J=1,NSCAT),(SCALE(I),I=1,NSF),IREJ,F
     &OMAX,FOMIN,RHOMN,RHOMX,(JSCRF(I),I=1,NSF),NRSF,WINKEL,AG,BG,CG    
     &,NPAR,(PAR(I),(PARNAM(J,I),J=1,2),IPAR(I),IRPAR(I),I=1,NPAR),FISGN
     &,IFSGN,EGGY,KEXTRF,BDF,GAMMAP,WVLNTH,INTPOL,JUNK6,INEUT,NWT,IWT   
     &,(INTER(I),I=1,NATOM),((Y(J,I),I=1,INTPOL),J=1,35),((Z(J,I),I=1,IN
     &TPOL),J=1,35),IUB,EXTRAN,IZACH,IREJEC,UNTEN,OBEN,MODELL,(ALF(I),I=
     &1,NATOM)                                                          
      READ(ISC1 )                                                       
     &((ANOM(J,I),J=1,2),I=1,NSCAT),IEXT,EXTIN,((TAG(J,I),J=1,2),ISOT(I)
     &,(ISCAT(J,I),J=1,2),OCCA(I),OCCUP(I),SITE(I),(XYZ(J,I),J=1,3),    
     &(BETA(J,I),J=1,6),(CR(J,I),J=1,10),(DR(J,I),J=1,15),I=1,NATOM),   
     &LDEGN,(KDEGN(I),I=1,LDEGN),HELPE,HELPF,IHIGH,((ER(J,I),J=1,21),(FR
     &(J,I),J=1,28),I=1,IHIGH) 
  131 read(isc1) itemp
      if(itemp) 132,131,132
  132 read(isc1) mdegn,(cnstr(i),i=1,mdegn)
      if(mdegn.gt.1000) type 133
  133 format(' TOO many constraints - increase array CNSTR !')
      if(mdegn.gt.1000) call exit
      j = mdegn + 1
      nmxl = (mrank*(mrank+1))/2 + mdegn
      read(isc1) (cv(i),i=j,nmxl)
      read(isc1) array(1)
      read(isc1) (bffd(i),i=1,mrank),(jirl(i),i=1,84*natom)
c      type 9133, (i,bffd(i),i=1,mrank)
c 9133 format(1h ,(10(i4,f8.4)))
c      type 9132, (jirl(i),i=1,84*natom)
c 9132 format(1h ,40(i3))
c
c---- calculate correlation matrix
c
c     ij = mdegn + 1
c     do 135 i=1,mrank
c     do 134 j=i,mrank
c     corr(i,j)=cv(ij)/(bffd(i)*bffd(j))
c 134 ij = ij + 1
c     type 9136, i,(j,corr(i,j),j=i,mrank)
c9136 format(5h0 ROW,i4/(10(i4,f8.4)))
c 135 continue
c
c---- save old parameters
c
      do 137 i =  1,natom
      do 1371 j = 1,3
 1371 xyzo(j,i) = xyz(j,i)
      do 1372 j = 1,6
 1372 betao(j,i) = beta(j,i)
      do 1373 j = 1,10
 1373 cro(j,i) = cr(j,i)
      do 1374 j = 1,15
 1374 dro(j,i) = dr(j,i)
      if(i.gt.6) goto 137
      do 1375 j = 1,21
 1375 ero(j,i) = er(j,i)
      do 1376 j = 1,28
 1376 fro(j,i) = fr(j,i)
  137 continue                          
C                                                                       
C---- PRINT COVER OF THE PROGRAM                                        
C                                                                       
      IF(IMARK.NE.0) GOTO 1112
      CALL SPACE(1,9,32,IOUT,NPAGE)                                     
      WRITE(IOUT,1111)LINE                                              
 1111 FORMAT(         28X,42H P R O M E T H E U S  HAS BEEN CALLED FOR  
     &  42H CALCULATION OF PROBABILITY DENSITY MAPS        /////            
     &40X,51H PROGRAMME FOR CALCULATION OF THE PROBABILITY DEN -  /     
     &40X,51H SITY FUNCTION OF THE EDGEWORTH, THE  GRAM-CHARLIER  /     
     &40X,51H AND THE SO CALLED ALPHA-FORMALISM  ( ALPHA=0 MEANS  /     
     &40X,51H GRAM - CHARLIER  SERIES EXPANSION ).  IN THE FIRST  /     
     &40X,51H CASE ALTERNATIVELY EXTENDED OR NORMAL  MAPS UP  TO  /     
     &40X,51H FOURTH ORDER CUMULANT PRODUCTS CAN BE SELECTED. IN  /     
     &40X,51H THE SECOND CASE NORMAL PDF-MAPS UP TO  SIXTH ORDER  /     
     &40X,51H CUMULANTS ARE IMPLEMENTED. INSTEAD OF THE PROBABI-  /     
     &40X,51H LITY  DENSITY THE  POTENTIAL  CAN  BE  CALCULATED.  /     
     &40X,51H INFORMATION  ABOUT  THIS  PROGRAMME  IS  AVAILABLE  /     
     &40X, 7H FROM  ,11A4)                                              
      WRITE(IOUT,240)LABEL                                              
  240 FORMAT(6(/),39X,41H FILE-LABEL OF THE LOGICAL UNIT MAIN =   ,     
     &3A4,10(/),35X,11HSUBSTANCE =,13X,13HTEMPERATURE =,13X, 9HSPECIAL =
     &)                                                                 
      CALL SPACE(2,9,32,IOUT,NPAGE)                                     
C                                                                       
C---- IS FOR THE USED STRUCTURE FACTOR EQUATION A PDF PROGRAMMED?       
C                                                                       
 1112 IF (MODELL.EQ.3) GOTO 912                                         
C                                                                       
C---- SAVE LENGTH OF CELL FOR OUTPUT                                    
C                                                                       
      CELLA=AG                                                          
      CELLB=BG                                                          
      CELLC=CG                                                          
C                                                                       
C---- TRANSFORM LATTICE PARAMETERS IN CORRESPONDENCE TO THE SELECTED    
C     SECTION.                                                          
C                                                                       
      CALL ARB(ISECT,CELL,ALFA,BET,WINKEL(3),AG,BG,CG,JACHS)            
C                                                                       
C---- SAVE GRID SELECTION FOR OUTPUT.                                   
C                                                                       
      MX  = IX                                                          
      MXI = IXI                                                         
      MXF = IXF                                                         
      MY  = IY                                                          
      MYI = IYI                                                         
      MYF = IYF                                                         
      MZ  = IZ                                                          
      MZI = IZI                                                         
      MZF = IZF                                                         
C                                                                       
C---- TRANSFORM INPUT PARAMETERS IN CORRESPONDENCE TO SELEC. SECTION.   
C                                                                       
      CALL ACHSE(ISECT,IX,IXI,IXF,IY,IYI,IYF,IZ,IZI,IZF)                
C                                                                       
C---- EXTENDED OR NORMAL MAP TO BE CALCULATED?                          
C                                                                       
      IPDF=IPDFF(6)                                                     
      IF(IPDF.EQ.0) GOTO 241                                            
C                                                                       
C---- EXTENDED REQUIRED. IS IT AN EDGEWORTH MAP?                        
C                                                                       
      IF(MODELL.NE.1) GOTO 908                                          
C                                                                       
C---- YES IT IS. READ FILE CONTAINING HARMONIC PARAMETERS.              
C                                                                       
      REWIND IHARMO                                                     
      READ  (IHARMO) ((XH(I,J),I=1,3),(BH(L,J),L=1,6),J=1,NATOM)        
      IEXT = 1                                                          
C                                                                       
C---- IS FIFTH AND/OR SIXTH ORDER CUMULANT COEF. DESIRED BUT NOT        
C     AVAILABLE?                                                        
C                                                                       
  241 IF(MODELL.EQ.2) GOTO 242                                          
      IF(IPDFF(4)+IPDFF(5)) 242,242,917                                 
C                                                                       
C---- PREPARE PARAMETERS FOR POINTARRAY OUTPUT, IF DESIRED.             
C                                                                       
  242 IF (NPNT.EQ.0) CALL SHERE                                         
C                                                                       
C---- PREPARE FILE S2, IF DATAS SHOULD SAVED FOR LATER PLOTTING.        
C     IN CASE OF SPECIAL PLANE RECORD WILL BEEN WRITTEN LATER.          
C                                                                       
      IF(IPLT.EQ.0 .OR. IPLANE.EQ.1) GOTO 101                           
      REWIND S2                                                         
      WRITE(S2)LABEL,TITLE,CELL,JACHS                                   
      REWIND S3                                                         
      WRITE(S3)LABEL,TITLE,CELL,JACHS                                   
      REWIND S4                                                         
      WRITE(S4)LABEL,TITLE,CELL,JACHS                                   
      REWIND S5                                                         
      WRITE(S5)LABEL,TITLE,CELL,JACHS                                   
      REWIND S6                                                         
      WRITE(S6)LABEL,TITLE,CELL,JACHS                                   
C                                                                       
C---- PRINT LABEL,TITLE AND STARTING TIME.                              
C                                                                       
  101 WRITE(IOUT,50)(TITLE(L),L=1,17),LABEL                             
   50 FORMAT(1H1,25HTITLE OF THE THIS PROGRAM /,1X,25(1H-)/1X,17A4//    
     &  35H0TITLE OF THE LEAST-SQUARES-PROGRAM /1X,34(1H-)/1X,3A4/)    
      WRITE(IOUT,43)                                                    
   43 FORMAT(14H0STARTING TIME/1X,13(1H-))                            
      CALL ANFANG(0)                                                    
C                                                                       
C---- PRINT PARAMETERS READ IN.                                         
C                                                                       
      WRITE(IOUT,41)                                                    
   41  FORMAT(//17H0INPUT-PARAMETERS/1X,16(1H-)/)                        
C                                                                       
C---- WHICH CUMULANTS SHOULD BE INCLUDED?                               
C                                                                       
      II=1                                                              
      DO 4111 I = 1,5                                                   
      IAUSG(I ) = JNOT                                                  
      IF(IPDFF(I).EQ.0) GO TO 4111                                      
      IAUSG(II) = JCUM(I)                                               
      II = II+1                                                         
 4111 CONTINUE                                                          
C                                                                       
C---- IS AT LEAST ONE CUMULANT SELECTED?                                
C                                                                       
      II = II-1                                                         
      IF(II.EQ.0) GOTO 612                                              
C                                                                       
C---- YES IT IS. ASSIGN TEXT CORRESPONDING TO IPUT PARAMETERS.          
C                                                                       
C---- SECTION TO BE SAVED FOR LATER PLOTTING?                           
  414 IDATA(1) = IDAT(IPLT+1)                                           
C---- POINT ARRAY TO BE PRINTED?                                        
      IDATA(2) = IDAT(NPNT+2)                                           
      IF(NPNT.NE.0) IDATA(2) = IDAT(NPNT)                               
C---- DOUBLE OR SINGLE SPACED OUTPUT?                                   
      IDATA(3) = IDAT(LSPC+3)                                           
      IDATA(4) = IDAT(5)                                                
C---- WHAT KIND OF STRUCTURE FACTOR EQUATION HAD BBEN USED?             
      DO 51 I=1,5                                                       
   51 IDATA(I+4)=ISER((MODELL-1)*5+I)                                   
C                                                                       
C---- IF SPECIAL PLANE SELECTED, CALCULATE PARAMETERS OF PLANE.         
C                                                                       
      IF(IPLANE.EQ.0) GOTO 511                                          
       CALL PLANE        (E1,E2,E3,MX,MY,MXI,MXF,MYI,MYF,ITYP1,         
     & G,XORG,YORG,ZORG,CELLA,CELLB,VECTR1,VECTR2,IERR)                 
C                                                                       
C---- WRITE FIRST RECORD FOR SPECIAL PLANE ON PLOTTING FILE.            
C                                                                       
      IF(IPLT.EQ.0) GOTO 511                                            
      ZERO = 0.                                                         
      ANGLE = 90.                                                       
      REWIND S2                                                         
      WRITE(S2)LABEL,TITLE,CELLA,CELLB,ZERO,ZERO,ZERO,ANGLE,JACHS       
      REWIND S3                                                         
      WRITE(S3)LABEL,TITLE,CELLA,CELLB,ZERO,ZERO,ZERO,ANGLE,JACHS       
      REWIND S4                                                         
      WRITE(S4)LABEL,TITLE,CELLA,CELLB,ZERO,ZERO,ZERO,ANGLE,JACHS       
      REWIND S5                                                         
      WRITE(S5)LABEL,TITLE,CELLA,CELLB,ZERO,ZERO,ZERO,ANGLE,JACHS       
      REWIND S6                                                         
      WRITE(S6)LABEL,TITLE,CELLA,CELLB,ZERO,ZERO,ZERO,ANGLE,JACHS       
C                                                                       
C---- HOW MANY ANGSTOEM PER GRID IN EACH DIRECTION.                     
C                                                                       
  511 IF(MX.EQ.0) MX = 1                                                
      SCX=CELLA/FLOAT(MX)                                               
      IF(MY.EQ.0) MY = 1                                                
      SCY=CELLB/FLOAT(MY)                                               
      IF(MZ.EQ.0) MZ = 1                                                
      SCZ=CELLC/FLOAT(MZ)                                               
      IF(IPLANE.EQ.1) MZ=0                                              
C                                                                       
C---- PRINT IT NOW.                                                     
C                                                                       
 512  WRITE(IOUT,52) (IDATA(4+I),I=1,5),IDAT(1 +IEXT),IDAT(IPOTET+1),   
     & ISEC(ISECT+1),(IDATA(I),I=1,4), GAM, EXXAG, IMINI, IMAXI,IDAT(   
     & IPOTET+1),IDAT(IPLANE+1),mc,sifa,MXI,MXF,MX,IORIGN(IOR+1),SCX,
     & MYI,MYF,MY,SCY,MZI,MZF,MZ,SCZ                                              
      mc = mc + 1
C                                                                       
   52 FORMAT           (17H0SERIES-EXPANSION,11X,8HEXTENDED,12X,9HPOTENT
     &IAL,12X,7HSECTION,13X,9HPLOT-FILE,15X,11HPOINT-ARRAY,//,1X,5A4,10X
     &,A4,16X,A4,16X,A4,17X,A4,20X,A4//                                 
     &11H0LINE-SPACE,17X,12HOUTPUT-ANGLE, 8X,12HEXXAGERATION, 9X,16HOMIT
     &TED - POINTS,4X,11HTEMPERATURE,13X,13HSPECIAL PLANE// 3X,2A4,18X, 
     &F8.2,12X,F10.3, 9X,I4,5H TILL,I5,10X,A4,21X,A4//                  
     &   14h0MC-ITERATIONS,4x,11hERROR-LEVEL,//,4x,i6,10x,f6.4,6h SIGMA,
     &//                    13H0GRID-SPACING,4X,3H>>>,8X,4HAXIS,16X,4HFR
     &OM,12X,4HTILL,12X,4HUNIT,12X,6HORIGIN,10X,15HANGSTROEMS/GRID//28X,
     &6HX-AXIS,9X,I9,12X,I4,12X,I4,13X,A4,11X,F9.4/28X,6HY-AXIS,9X,I9,12
     &X,I4,12X,I4,28X,F9.4,/,28X,6HZ-AXIS,9X,I9,12X,I4,12X,I4,28X,F9.4//
     &//)                                                               
C                                                                       
      WRITE(IOUT,525)IAUSG                                              
 525  FORMAT(32H0INCLUDED CUMULANTS ARE OF RANK ,15X,5(A7,4X)///)       
C                                                                       
C---- IN CASE OF SPECIAL PLANE OUTPUT WILL BE DIFFERENT FROM NORM. PLANE
C                                                                       
      IF(IPLANE) 5252,5251,5252                                         
C                                                                       
C---- NORMAL PLANE.                                                     
 5251 WRITE(IOUT,53)                                                    
      GOTO 5254                                                         
C---- SPECIAL PLANE.                                                    
 5252 WRITE(IOUT,531)                                                   
      J=0                                                               
C---- MORE THAN ONE ATOM SELECTED?                                      
      DO 5253 I = 1,NATOM                                               
 5253 IF(IPLOT(I).NE.0) J = J + 1                                       
      IWRITE=0                                                          
      IF(J.EQ.1) IWRITE = 1                                             
C                                                                       
   53 FORMAT(30H0PDF-CALCULATION FOR THE ATOMS,3X,3H>>>,4X,5HLABEL,15X, 
     &8HPOSITION,20X,16HGRID-COORDINATES, 8X,13HT-FACTOR TYPE//)        
  531 FORMAT(30H0PDF-CALCULATION FOR THE ATOMS,3X,3H>>>,4X,5HLABEL,15X, 
     &8HPOSITION,12X,23HXGRID  YGRID   DISTANCE, 7X,13HT-FACTOR TYPE//) 
C                                                                       
 5254 DO 54 I=1,NATOM                                                   
C                                                                       
C---- ARE THERE ATOMS BEEN SELECTED?                                    
C                                                                       
      IF(IPLOT(I).EQ.0)GOTO 54                                          
C                                                                       
C---- YES. CHECK FOR LEGAL SELECTION OF CUMULANTS.                      
C                                                                       
      IF(ISOT(I).LE.1 .AND. IPDFF(2).NE.0) IPLOT(I)=0                   
      IF(ISOT(I).LE.2 .AND. IPDFF(3).NE.0) IPLOT(I)=0                   
      IF(ISOT(I).LE.3 .AND. IPDFF(4).NE.0) IPLOT(I)=0                   
      IF(ISOT(I).LE.4 .AND. IPDFF(5).NE.0) IPLOT(I)=0                   
C                                                                       
      IF(IPLANE) 5256,5255,5256                                         
C                                                                       
C---- NORMAL PLANE                                                      
C                                                                       
 5255 F1=FLOAT(IOR)                                                     
C                                                                       
C---- PRINT SELECTED ATOMS.                                             
C                                                                       
      XH1 = XH(1,I)*IX*F1
      XH2 = XH(2,I)*IY*F1
      XH3 = XH(3,I)*IZ*F1
      XYZ1 = XYZ(1,I)*IX*F1
      XYZ2 = XYZ(2,I)*IY*F1
      XYZ3 = XYZ(3,I)*IZ*F1
      ITP1 = ITYPE(ISOT(I)*2+1)
      ITP2 = ITYPE(ISOT(I)*2+2)
      IF(IEXT.EQ.1)WRITE(IOUT,55)(TAG(J,I),J=1,2),(XH(J,I),J=1,3),  
     &XH1,XH2,XH3,IPT1,IPT2
      IF(IEXT.EQ.0)                                                     
     & WRITE(IOUT,55)(TAG(J,I),J=1,2),(XYZ(J,I),J=1,3),XYZ1,  
     &XYZ2,XYZ3,ITP1,ITP2
      GOTO 54                                                           
C                                                                       
C---- SPECIAL PLANE.                                                    
C                                                                       
 5256 IF(IEXT) 5257,5258,5257                                           
C                                                                       
C---- EXTENDED MAP. (IF IERR=1 GOTO 9003 TO WRITE ERRORMESSAGE AND EXIT)
 5257 IF(IERR.EQ.1) GOTO 9003                                           
      XHELP(1)=XH(1,I)                                                  
      XHELP(2)=XH(2,I)                                                  
      XHELP(3)=XH(3,I)                                                  
      CALL ATOPLA(E1,E2,E3,XHELP,MX,MY,G,POUT)                          
      WRITE(IOUT,5259)(TAG(J,I),J=1,2),(XH(J,I),J=1,3),POUT,(ITYPE(     
     & ISOT(I)*2+L),L=1,2)                                              
 5259 FORMAT(42X,2A3,9X,3(F5.3,2X),2X,2(F5.1,2X),F6.2,8X,2A4)           
      GOTO 54                                                           
C                                                                       
C---- NORMAL MAP (SPECIAL PLANE).                                       
 5258 IF(IERR.EQ.1) GOTO 9003                                           
      XHELP(1)=XYZ(1,I)                                                 
      XHELP(2)=XYZ(2,I)                                                 
      XHELP(3)=XYZ(3,I)                                                 
      CALL ATOPLA(E1,E2,E3,XHELP,MX,MY,G,POUT)                          
      WRITE(IOUT,5259)(TAG(J,I),J=1,2),(XYZ(J,I),J=1,3),POUT,(ITYPE(    
     & ISOT(I)*2+L),L=1,2)                                              
   54 CONTINUE                                                          
C                                                                       
   55 FORMAT(42X,2A3,9X,3(F5.3,2X),8X,3(F5.1,2X), 8X,2A4)               
  120 CONTINUE                                                          
C                                                                       
C---- CORRESPONDING TO THE NON ZERO VALUES OF IPDFF CUMULANTS OF        
C     ANY RANK ARE INCLUDED.                                            
C     IPDFF(1) = 1 MEANS INCLUSION OF TENOR OF    RANK TWO.             
C     IPDFF(2) = 2 MEANS INCLUSION OF TENSOR OF RANK THREE.             
C     AND SO ON.                                                        
C                                                                       
      ISKEW   =  1.                                                     
      ITHREE  =  1                                                      
      IFOUR   =  1                                                      
      IFIVE   =  1                                                      
      ISIX    =  1                                                      
      ICUM1   =  1                                                      
      ICUM2   =  1                                                      
C                                                                       
      IF(IPDFF(1).NE.0) GOTO 13                                         
C---- DO NOT INCLUDE THE SECOND RANK CUMULANT.                          
      ISKEW   = 0.                                                      
   13 IF(IPDFF(2).NE.0) GOTO 14                                         
C---- DO NOT INCLUDE THE THIRD RANK TENSOR.                             
      ICUM1   = 0                                                       
      ITHREE  = 0                                                       
   14 IF(IPDFF(3).NE.0) GOTO 15                                         
C---- DO NOT INCLUDE THE FOURTH RANK TENSOR.                            
      ICUM2   = 0                                                       
      IFOUR   = 0                                                       
   15 IF(IPDFF(4).NE.0) GOTO 16                                         
C---- DO NOT INCLUDE THE FIFTH RANK TENSOR.                             
      IFIVE   = 0                                                       
   16 IF(IPDFF(5).NE.0) GOTO 17                                         
      ISIX    = 0                                                       
   17 CONTINUE                                                          
C                                                                       
C---- IN CASE OF POTENTIAL CALCULATION A SELECTION OF CUMULANTS         
C     WITHOUT  SECOND ORDER IS MEANINGLESS.                             
C                                                                       
      IF(IPOTET.EQ.0 .OR. ISKEW .GT. 0.) GOTO 18                        
      ILINE=12                                                          
      WRITE(IOUT,923)                                                   
      ISKEW=1.                                                          
C                                                                       
C---- CALCULATE GRID SPACING.                                           
C                                                                       
  18  IF(IPLANE)181,181,182                                             
 181  XX=1./FLOAT(IX)                                                   
      YY=1./FLOAT(IY)                                                   
      ZZ=1./FLOAT(IZ)                                                   
 182  IXA=IXF-IXI+1                                                     
      IYA=IYF-IYI+1                                                     
      IZA=IZF-IZI+1                                                     
C                                                                       
C---- SHOULD THE OUTPUT ARRAY SINGLE OR DOUBLE SPACED?                  
C                                                                       
      LSPCH=LSPC                                                        
      IF(LSPC.NE.0) GOTO 2222                                           
      LSPC=ISPC                                                         
      GOTO 3333                                                         
 2222 LSPC=IBNK                                                         
 3333 CONTINUE                                                          
C                                                                       
C---- CALCULATION OF NECESSARY AND AVAILABLE CORE ALLOCATIONS.          
C                                                                       
      ICORE=IXA*IYA                                                     
      WRITE(IOUT,4444)ICORE,MAXIMA                                      
 4444 FORMAT(17H0REQUIRED STORAGE  ,/,1X,16(1H-)//,1X,   43HUSED ELEMENT
     &S OF THE ARRAY NAMED "ARRAY"    ,I6,1H.,15X,9HAVAILABLE,I6,1H.)   
      IF(ICORE.GT.MAXIMA) GOTO 910                                      
      ILINE=49+NATOM+ILINE                                              
C                                                                       
C---- WARNING IN CASE OF EDGEWORTH CALCULATIONS.                        
C                                                                       
      IF(ILINE.LT.58 .AND. 58 .LT. NPAGE) GO TO 9001                    
      ILINE=1                                                           
      WRITE(IOUT,9002)                                                  
 9002 FORMAT(1H1)                                                       
 9001 IF(MODELL.EQ.1)WRITE(IOUT,19)                                     
   19 FORMAT(//18H0IMPORTANT NOTICE    /1X,16(1H-)//                    
     &61H0THE EDGEWORTH-SERIES EXPANSION USED IN THIS PROGRAM IS A BAD /
     &61H0BUT THE ONLY KNOWN APPROXIMATION FOR THE FOURIER - TRANSFORM /
     &61H0OF THE TEMPERATUREFACTOR. (MODEL NUMBER 1)                   /
     &61H0BE CAREFUL IN THE INTERPRETATION OF THE SKEW-MAP. THE EDGE - /
     &61H0WORTH-SERIES IS IN THIS CASE   N O T  DIFFERENT TO THE GRAM- /
     &27H0CHARLIER-SERIES EXPANSION. )                                  
C                                                                       
C---- WRITE INFO. ABOUT SPECIAL PLANE IF DESIRED.                       
C                                                                       
 9003 IF(IPLANE.NE.0) CALL SRBPLA(G, E1,E2,E3, XORG,YORG,ZORG, VECTR1,  
     & VECTR2,IOUT,IERR,ITYP1,NPAGE,POUT,MX,MY,IWRITE)                  
C                                                                       
C---- WRITE LIST OF TENSOR-ELEMENTS OF EACH RANK.                       
C                                                                       
      JUNK=0                                                            
      do 5 i=1,10000
      qp(i) = 0.0
      qm(i) = 0.0
      iplus(i) = 0
      iminus(i) = 0
    5 continue
      jpri = 0
C                                                                       
c
c---- s t a r t    B I G  M o n t e - C a r l o  L o o p ----
c
      do 2500 imoca = 1,mc
      imc = imoca - 1
c---- initialize counter for special plane
      jndex = 0
      type 301, imc
  301 format(' imc = ',i5)
c
      DO 300 I=1,NATOM                                                  
      IF(IPLOT(I).LE.0) GOTO 300                                        
      JUNK=1                                                            
C                                                                       
C---- IN CASE OF ISOTROPIC T-FACTOR CHANGE TO ANISOTROPIC FORM.         
C                                                                       
      IF(ISOT(I).NE.0) GOTO 319                                         
      TEMP=BETA(1,I)                                                    
      BETA(6,I)=0.25*RECELL(2)*RECELL(3)*RECELL(4)*TEMP                 
      BETA(5,I)=0.25*RECELL(1)*RECELL(3)*RECELL(5)*TEMP                 
      BETA(4,I)=0.25*RECELL(1)*RECELL(2)*RECELL(6)*TEMP                 
      BETA(3,I)=0.25*RECELL(3)**2*TEMP                                  
      BETA(2,I)=0.25*RECELL(2)**2*TEMP                                  
      BETA(1,I)=0.25*RECELL(1)**2*TEMP                                  
C                                                                       
  319 if(imc.ne.0) goto 300
      IF(MODELL.NE.1) GOTO 500                                          
C                                                                       
C---- CALCULATION OF AN EDGEWORTH SERIES EXPANSION.                     
C     (THE EXISTENCE OF TWO OUTPUT PARTS IS HISTORICALLY)               
C                                                                       
       WRITE(IOUT,401)                                                  
C                                                                       
       DO 310 JK=1,LOOP(ISOT(I)+1)                                      
          DO 320 JL=1,8                                                 
 320      IOUTP(JL)=IBNK                                                
       IF(JK.GT.1)GOTO 330                                              
       OUTP(1)=TAG(1,I)                                                 
       OUTP(2)=TAG(2,I)                                                 
 330   IF(JK.GT.3) GOTO 340                                             
       IF(IEXT.EQ.1) IOUTP(3)=IPRI(JK)                                  
                     IOUTP(4)=IPRI(JK)                                  
 340   IF(JK.GT.6) GOTO 350                                             
       IF(IEXT.EQ.1) IOUTP(5)=IPRI(3+JK)                                
                     IOUTP(6)=IPRI(3+JK)                                
 350   IF(JK.GT.10) GOTO 360                                            
                     IOUTP(7)=IPRI(9+JK)                                
 360   IF(ISOT(I).LT.3) GOTO 370                                        
                     IOUTP(8)=IPRI(19+JK)                               
 370   IF(JK.GT.3) GOTO 380                                             
       IF(IEXT.EQ.1)WRITE(IOUT,403)IOUTP(3),XH(JK,I),IOUTP(5),BH(JK,I)  
       IF(IEXT.EQ.0)WRITE(IOUT,412)                                     
       WRITE(IOUT,404)(OUTP(L),L=1,2),IOUTP(4),XYZ(JK,I),IOUTP(6),      
     & BETA(JK,I),IOUTP(7),CR(JK,I)                                    
       IF(ISOT(I).LT.3) GOTO 380                                        
       WRITE(IOUT,405)IOUTP(8),DR(JK,I)                                 
       GOTO 310                                                         
 380   IF(JK.GT.6)GOTO 390                                              
       IF(IEXT.EQ.1) WRITE(IOUT,406)IOUTP(5),BH(JK,I)                   
       IF(IEXT.EQ.0) WRITE(IOUT,412)                                    
       WRITE(IOUT,407)IOUTP(6),BETA(JK,I),IOUTP(7),CR(JK,I)             
       IF(ISOT(I).LT.3) GOTO 310                                        
       WRITE(IOUT,405) IOUTP(8),DR(JK,I)                                
       GOTO 310                                                         
 390   IF(JK.GT.10) GOTO 400                                            
       WRITE(IOUT,408) IOUTP(7),CR(JK,I)                                
       IF(ISOT(I).LT.3) GOTO 310                                        
       WRITE(IOUT,405) IOUTP(8),DR(JK,I)                                
       GOTO 310                                                         
 400   IF(ISOT(I).LT.3) GOTO 310                                        
       WRITE(IOUT,409) IOUTP(8),DR(JK,I)                                
  310  CONTINUE                                                         
       GOTO 300                                                         
C                                                                       
C----  OUTPUT FOR GRAM CHARLIER EXPANSION AND ALPHA FORM.(GC)           
C                                                                       
 500   CONTINUE                                                         
C                                                                       
C----  PRINT HEADING.                                                   
C                                                                       
       IF(MODELL.EQ.2) WRITE(IOUT,600)                                  
       IF(MODELL.EQ.4) WRITE(IOUT,608)                                  
C                                                                       
       DO 510 JK=1,LOOP(ISOT(I)+1)                                      
C                                                                       
          WRITE(IOUT,412)                                               
          IF(JK.GT.21) GOTO 585                                         
          IF(JK.GT.15) GOTO 584                                         
          IF(JK.GT.10) GOTO 583                                         
          IF(JK.GT. 6) GOTO 582                                         
          IF(JK.GT. 3) GOTO 581                                         
          IF(JK.GT. 1) GOTO 580                                         
C                                                                       
C---- PRINT LABEL                                                       
C                                                                       
          IF(MODELL.EQ.2) WRITE(IOUT,601)TAG(1,I),TAG(2,I)              
          IF(MODELL.EQ.4) WRITE(IOUT,609)TAG(1,I),TAG(2,I),ALF(I)       
C                                                                       
C---- PRINT X-COORDINATES.                                              
C                                                                       
  580     WRITE(IOUT,602) IPRI(JK),XYZ(JK,I)                            
C                                                                       
C---- PRINT B-TENSOR, C-TENSOR, D-TENSOR AND SO ON.                     
C                                                                       
  581 BETA(JK,I) = BETA(JK,I)*100.
          WRITE(IOUT,603) IPRI(JK+ 3),BETA(JK,I)                   
      BETA(JK,I) = BETA(JK,I)/100.
  582 IF(ISOT(I).GT.1)                                                  
     &    WRITE(IOUT,604) IPRI(JK+ 9),  CR(JK,I)                        
  583 IF(ISOT(I).GT.2)                                                  
     &    WRITE(IOUT,605) IPRI(JK+19),  DR(JK,I)                        
  584 IF(ISOT(I).GT.3)                                                  
     &    WRITE(IOUT,606) IFUENF (JK),  ER(JK,I)                        
  585 IF(ISOT(I).GT.4)                                                  
     &    WRITE(IOUT,607) ISECHS (JK),  FR(JK,I)                        
C                                                                       
  510  CONTINUE                                                         
C                                                                       
 300  CONTINUE                                                          
C                                                                       
C---- FORMAT STATEMENTS FOR PARAMETER OUTPUT LIST.                      
C                                                                       
 401  FORMAT(1H1,5HLABEL,9X,11HX(HARMONIC),10X,9HX -TENSOR,10X,         
     & 11HB(HARMONIC),10X,9HB -TENSOR, 6X,15HC-TENSOR(*1000),4X,        
     & 16HD-TENSOR(*10**4) //)                                          
 403  FORMAT(1H0,13X,A2,F10.5,27X,A4,F9.5)                              
 404  FORMAT(1H+,2A4,26X,A2,F9.5,29X,A4,F7.5,6X,A4,1H=,F9.5)            
 405  FORMAT(1H+,108X,1HD,A4,1H=,F9.5)                                  
 406  FORMAT(1H0,52X,A4,F9.5)                                           
 407  FORMAT(1H+,74X,A4,F7.5,6X,A4,1H=,F9.5)                            
 408  FORMAT(1H0,91X,A4,1H=,F9.5)                                       
 409  FORMAT(1H0,108X,1HD,A4,1H=,F9.5)                                  
 412  FORMAT(1H0)                                                       
C                                                                       
 600  FORMAT(1H1,5HLABEL,     15X,10HX - COORD.,10X,10HB - TENSOR,10X,  
     &10HC - TENSOR,10X,10HD - TENSOR,10X,10HE - TENSOR,10X,7HF - TEN,  
     &3HSOR,   /,42X, 8H(*10**2)  ,12X, 8H(*10**3)  ,12X, 8H(*10**4)  , 
     &           12X, 8H(*10**5)  ,12X, 8H(*10**6)  ,/)                 
 601  FORMAT(1H+,2A4)                                                   
 602  FORMAT(1H+,20X,A2,F8.5)                                           
 603  FORMAT(1H+,    38X    ,A4,F8.5)                                   
 604  FORMAT(1H+,          58X      ,A4,1H=,F8.5)                       
 605  FORMAT(1H+,               77X             ,1HD,A4,1H=,F8.5)       
 606  FORMAT(1H+,                   96X         ,A7,F8.5        )       
 607  FORMAT(1H+,                         115X  ,A7,1H=,F8.5    )       
 608  FORMAT(1H1,5HLABEL,     15X,10HX - COORD.,10X,10HB - TENSOR,10X,  
     &10HC - TENSOR,10X,10HD - TENSOR,10X,10H  ALPHA   /                
     & 42X, 8H(*10**2)  ,12X, 8H(*10**3)  ,12X, 8H(*10**4)/)            
 609  FORMAT(1H+,2A4,105X,F6.3)                                         
C                                                                       
C---- ERRORMESSAGE, IN TWO CASES.  FIRST= NO ATOM FOR PDF-CALCULATION   
C                                         HAD BEEN SELECTED.            
C                                  SECOND = A PDF TYPE HAVE BEEN SELEC- 
C                                         TED, WHICH IS NOT AVAILABLE.  
C                                                                       
      IF(JUNK.EQ.1) GOTO 410                                            
      WRITE(IOUT,411)                                                   
  411 FORMAT(////13H0ERRORMESSAGE/1X,12(1H-)//                          
     &77H0ONLY ATOMS WITH INCONSISTENT TEMPERATURE FACTORS OR NO ATOM AT
     & ALL SELECTED. /31H I'M SORRY, I 'VE TO CALL EXIT.)               
      CALL EXIT                                                         
C                                                                       
 410  CONTINUE                                                          
C                                                                       
C---- LOOP THROUGH SELECTED ATOMS.                                      
C                                                                       
      DO 2400 IATOM =1,NATOM                                            
      IF(IPLOT(IATOM).LE.0) GOTO 2400
c
c---- some prelim calculation for mc
c
      if(imc.ne.0) goto 916
c
c---- calculate number of free Mc-variables
      mcfp = 0
      jirmin = 1000
      jirmax = 0
      do 914 i=1,84
      if(jir(i,iatom).eq.0) goto 914
      if(jirmin.gt.jir(i,iatom)) jirmin = jir(i,iatom)
      if(jirmax.lt.jir(i,iatom)) jirmax = jir(i,iatom)
      mcfp = mcfp + 1
  914 continue
c---- repak constraint information for use in RESET
      call repak(natom,mdegn)
      nij = (mcfp*(mcfp+1))/2
      call crack(jirmin,jirmax,mcfp,vcv,c,cl,mdegn,mrank)
c---- calculate KPARA (indexing variable for routine RESET)
      kpara = 0
      do 915 i=1,iatom-1
      jpara = kpara + 5
      if(isot(i).ge.1) jpara = jpara+5
      if(isot(i).ge.2) jpara = jpara+10
      if(isot(i).ge.3) jpara = jpara+15
      if(isot(i).ge.4) jpara = jpara+21
      if(isot(i).ge.5) jpara = jpara+28
  915 kpara = jpara
  916 IHOLD=IEXT                                                        
      IF(ISOT(IATOM).EQ.0) IEXT=0                                       
C                                                                       
C---- CHECK FOR ZERO CUMULANT PRODUCTS.                                 
C                                                                       
      NO3=0                                                             
      NO4=0                                                             
      NO5=0                                                             
      NO6=0                                                             
C                                                                       
      IAWAY=0                                                           
C                                                                       
      IF(ISOT(IATOM).EQ.1) GOTO 420                                     
C                                                                       
C---- IS THE THIRD CUMULANT PRODUCT NECESSARY FOR DESIRED PDF-MAP?      
C     IN CASE OF AN EDGEWORTH SERIES  FOURTH ORDER EQUATION  NEEDS      
C     TENSOR OF THIRD ORDER.                                            
C                                                                       
      IF(ITHREE.EQ.0 .AND. MODELL.NE.1) GOTO 4201                       
C                                                                       
C----  AND IS THIRD ORDER AVAILABLE?                                    
C                                                                       
      IF(ISOT(IATOM).LT.2) GOTO 420                                     
C                                                                       
C---- YES| IS IT A NULL PRODUCT?                                        
C                                                                       
      NO3=0                                                             
      SH=0.0                                                            
      DO 421 IE=1,10                                                    
  421 SH=SH+ABS(CR(IE,IATOM))                                           
      IF(SH.LT.1.E-11) GOTO 4200                                        
C---- NO, WHICH ELEMENTS ARE EXACTLY ZERO?                              
      NO3=1                                                             
      DO 4199 IE=1,10                                                   
      JTHREE(IE)=0                                                      
 4199 IF(ABS(CR(IE,IATOM)).GT.1.E-9) JTHREE(IE)=1                       
C                                                                       
C---- ONLY IN CASE OF EDGEWORTH SERIES.                                 
C     SHOULD THIRD ORDER TO BE OMITTED? (FOURTH ORDER NEEDS C(IJK))     
C                                                                       
 4200 IF(MODELL.NE.1) GOTO 4201                                         
      IF(NO3.EQ.1 .AND. ITHREE.EQ.0) IAWAY=1                            
C                                                                       
C---- HAS THE ATOM A FOURTH ORDER CUMULANT PRODUCT?                     
C                                                                       
 4201 IF(IFOUR.EQ.0) GOTO 4202                                          
C                                                                       
C---- IS FOURTH CUMULANT PRODUCT AVAILABLE.                             
C                                                                       
      IF(ISOT(IATOM).LT.3) GOTO 420                                     
C                                                                       
C---- IS THE FOURTH RANK TENSOR PRODUCT A NULL PRODUCT?                 
C                                                                       
      NO4=0                                                             
      SH=0.0                                                            
      DO 422 IE=1,15                                                    
  422 SH=SH+ABS(DR(IE,IATOM))                                           
      IF(SH.LT.1.E-11) GOTO 4202                                        
      NO4=1                                                             
C---- NO, WHICH ELEMENTS ARE EXACTLY ZERO?                              
      DO 4220 IE=1,15                                                   
      JFOUR(IE)=0                                                       
 4220 IF(ABS(DR(IE,IATOM)).GT.1.E-10)JFOUR(IE)=1                        
C                                                                       
C---- IS THE TEMPERATURE FACTOR OF AT LEAST FIFTH ORDER?                
C                                                                       
 4202 IF(ISOT(IATOM).EQ.3) GOTO 420                                     
C                                                                       
C---- IS FIFTH ORDER CUMULANT PRODUCT A NULL PRODUCT?                   
C                                                                       
      IF(IFIVE.EQ.0) GOTO 4203                                          
      NO5=0                                                             
      SH=0.0                                                            
      DO 423 IE=1,21                                                    
  423 SH=SH+ABS(ER(IE,IATOM))                                           
      IF(SH.LT.1.E-7) GOTO 4203                                         
      NO5=1                                                             
C                                                                       
C---- NOT A NULL PRODUCT. WHICH ELEMENTS OF THE TENSOR OF FIFTH RANK    
C     ARE EXACTLY ZERO.                                                 
C                                                                       
      DO 4631 IE=1,21                                                   
      JFIVE(IE)=0                                                       
 4631 IF(ABS(ER(IE,IATOM)).GT.1.E-7) JFIVE(IE)=1                        
C                                                                       
C---- IS THE TEMPERATURE FACTOR OF SIXTH ORDER?                         
C                                                                       
 4203 IF(ISOT(IATOM).EQ.4) GOTO 420                                     
C                                                                       
C---- IS SIXTH ORDER CUMULANT PRODUCT A NULL PRODUCT?                   
C                                                                       
      IF(ISIX.EQ.0) GOTO 420                                            
      NO6=0                                                             
      SH=0.0                                                            
      DO 4232 IE=1,28                                                   
 4232 SH=SH+ABS(FR(IE,IATOM))                                           
      IF(SH.LT.1.E-7) GOTO 420                                          
      NO6=1                                                             
C                                                                       
C---- NOT A NULL PRODUCT. WHICH ELEMENTS OF THE TENSOR OF SIXTH RANK    
C     ARE EXACTLY ZERO.                                                 
C                                                                       
      DO 4233 IE=1,28                                                   
      JSIX(IE)=0                                                        
 4233 IF(ABS(FR(IE,IATOM)).GT.1.E-7) JSIX(IE)=1                         
C                                                                       
C---- CALCULATE FLAG FOR JUMP OVER ANHARMONIC PART.                     
C                                                                       
 420  ITOT=NO3+NO4+NO5+NO6                                              
C                                                                       
C---- REMEMBER. IF YOU CALCULATE ONLY SECOND ORDER PDF, EXTENSION       
C     IS MEANINGLESS FOR COMPUTATION.                                   
C---- NOW CHECK THAT NOT A NULL RANK TENSOR SHOULD BE CALCULATED.       
C                                                                       
      IF(ITOT.EQ.0 .AND. (ITHREE+IFOUR+IFIVE+ISIX).NE.0)                
     & GOTO 921                                                         
C                                                                       
C---- CLEAR ORIGIN TRANSLATION VECTOR.                                  
C                                                                       
      DO 989 I=1,3                                                      
 989  R1(I)=0.0                                                         
C                                                                       
C---- SHUFFEL IATOM TO II TO PREVENT COMPILER WARNINGS.                 
      II=IATOM
c
c---- get random vector and correlated rv
c
      if(imc.eq.0) goto 4222
      call ranvec(mcfp,ranr,seed,twopi)
      call rancor(mcfp,nij,cl,ranr,ranc)
      call modify(mcfp,jirmin,jirmax,ii,ranc)
      call reset(ii,kpara,mdegn)
C                                                                       
C---- IF AN EDGEWORTH MAP TO BE CALCULATED CREATE TRANSFORMATIONS       
C     AND PARAMETERS.                                                   
C                                                                       
C---- INITIALIZE FLAG FOR   NO DIFFERENCE BETWEEN NORMAL AND EXTENDED   
C     MAP.                                                              
C                                                                       
 4222 IEMAP1=0                                                          
      IEMAP2=0                                                          
C                                                                       
C---- NO EXTENSION OR NO ANHARMONICITY?                                 
C                                                                       
      IF(IEXT.EQ.0 .OR.ITOT.EQ.0 ) GOTO 21                              
C                                                                       
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++   
C                                                                       
CEXTENDED MAP  EXTENDED MAP  EXTENDED MAP  EXTENDED MAP  EXTENDED MAP   
C                                                                       
C---- TRANSFORM TENSOR COMPONENTS CORRESPONDING TO SECTION SELECTION.   
C                                                                       
      IF(ISECT.EQ.0) GOTO 22                                            
      CALL TRANSF(ISECT,INDEX1,INDEX2,INDEX3,MIROR3)                    
C                                                                       
C---- DEVIATION OF TRANSLATION VECTOR DUE TO EXTENDED MAP.              
C                                                                       
   22 DO 23 I=1,3                                                       
   23 R1(I)=XYZ(I,II)-XH(I,II)                                          
C                                                                       
C---- CHANGE OF TENSOR OF RANK 2 DUE TO A EXTENDED MAP CALCULATION.     
C                                                                       
      IF(ISOT(II).LT.3 .OR. NO4.EQ.0) GOTO 2411                         
      DO 24 I=1,6                                                       
   24 R2(I)=(BETA(I,II)-BH(I,II))/TWOPI2                                
C                                                                       
C---- CALCULATE INVERSE MATRIX OF HARMONIC TEMPERATURE PART.            
C                                                                       
 2411 CALL UMKEHR(BH,II,QDET)                                           
C                                                                       
C                                                                       
C---- ARE THERE ANY CHANGES IN THE HARMONIC FACTOR DUE TO  EXTENSION OF 
C     MAP.                                                              
C                                                                       
      IF(ICUM1.EQ.0) GOTO 4231                                          
      IEMAP=  1                                                         
      SH=ABS(R1(1))+ABS(R1(2))+ABS(R1(3))                               
      IF(SH.LT.1.E-7) IEMAP1 = 0                                        
C                                                                       
C---- IS THERE ANY CHANGE ON THE SECOND RANK TENSOR DUE TO EXTENSION OF 
C     THE MAP AND IS EXTENDED MAP OF FOURTH ORDER?                      
C                                                                       
 4231 IF(ISOT(II).LT.3 .OR. ICUM2.EQ.0  ) GOTO 25                       
C                                                                       
      SH=0.0                                                            
      IEMAP2=1                                                          
      DO 426 I=1,6                                                      
  426 SH=SH+ABS(R2(I))                                                  
      IF(SH.LT.1.E-7) IEMAP2 = 0                                        
      GOTO 25                                                           
C                                                                       
C END OF FIRST PART OF EXTENDED MAP CALCULATION   END OF FIRST PART OF  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
C                                                                       
C                                                                       
C---- SAME FOR NORMAL MAP AND GRAM-CHARLIER SERIES EXPANSION -------    
C                                                                       
   21 IF(ISECT.EQ.0) GOTO 26                                            
      CALL TRANSF(ISECT,INDEX1,INDEX2,INDEX3,MIROR3)                    
C                                                                       
C---- MATRIX INVERSE OF HARMONIC PART OF TEMPERATURE FACTOR.            
C                                                                       
   26 CALL UMKEHR(BETA,II,QDET)                                         
C                                                                       
C---- CALCULATE SMALL FORMULAS OFTEN USED IN THE FOLLOWING.             
C     (NOT FOR MODEL NUMBER 4 (ALPHA FORMALISM (GC))                    
C                                                                       
   25 IF(MODELL.EQ.4) GOTO 251                                          
      CALL PREPAR(MODELL,II,CR)                                         
C                                                                       
C---- CALCULATE FLAGS FOR JUMPS                                         
C                                                                       
  251 JUMP1=1                                                           
      JUMP2=1                                                           
      JUMP3=1                                                           
      IF(NO3.EQ.0 .OR. IEMAP1.EQ.0 ) JUMP2=0                            
      IF(NO4.EQ.0 .OR. JUMP2.EQ.0  ) JUMP1=0                            
      IF(NO4.EQ.0 .OR. IEMAP1.EQ.0 .AND. IEMAP2.EQ.0 ) JUMP3=0          
C---- IN CASE OF AN EXTENDED EDGEWORTH MAP ALL HERMITE POLYNOMS OF 4 TH 
C     ORDER CAN BE NECESSARY. THEREFORE CALCULATE THESE.                
      IF(JUMP1*JUMP3)100,100,99                                         
   99 DO 102 IE=1,15                                                    
  102 JFOUR(IE)=1                                                       
  100 CONTINUE                                                          
C                                                                       
C---- ROMAXT SHOULD BE MAXIMUM OF CALCULATED PDF OF FIRST SECTION.      
      ROMAXT = - 1000000.                                               
C                                                                       
C POINT BY POINT CALCULATION OF THE PROBABILITY MAP IN DESIRED SECTIONS.
C                                                                       
C ========= LOOP OVER SECTIONS =========================================
C                                                                       
       DO 890  JZSEC=1,IZA                                              
       ROMIN = 100000.                                                  
       ROMAX =-100000.                                                  
       jdruck = jdr(1)
       IZAEHL= 0                                                        
C                                                                       
C----  IN CASE OF A SPECIAL PLANE CALC. X,Y,Z WILL BE CALC. BY CORD.    
       IF(IPLANE)891,891,892                                            
  891  X(INDEX3)=ZZ*FLOAT(JZSEC+IZI-1)                                  
       IF(IOR.NE.0) X(INDEX3)=X(INDEX3)-XYZ(INDEX3,II)-R1(INDEX3)       
       X(INDEX3) = MIROR3*X(INDEX3)                                     
  892  CONTINUE                                                         
C                                                                       
C   ========== LOOP OVER VERTICAL AXIS =================================
C                                                                       
        DO 2200 JYSEC=1,IYA                                             
C                                                                       
        IF(IPLANE) 2201,2201,2202                                       
 2201   X(INDEX2)=YY*FLOAT(JYSEC+IYI-1)                                 
        IF(IOR.NE.0) X(INDEX2)=X(INDEX2)-XYZ(INDEX2,II)-R1(INDEX2)      
 2202   CONTINUE                                                        
C                                                                       
C      ====== LOOP OVER HORIZONTAL AXIS ================================
C                                                                       
         DO 2100 JXSEC=1,IXA                                            
C                                                                       
         IF(IPLANE) 2101,2101,2102                                      
 2101    X(INDEX1)=XX*FLOAT(JXSEC+IXI-1)                                
         IF(IOR.NE.0) X(INDEX1)=X(INDEX1)-XYZ(INDEX1,II)-R1(INDEX1)     
         GOTO 2103                                                      
C                                                                       
C----- SPECIAL COORDINATES BY ROUTINE CORD.                             
C                                                                       
 2102    CALL CORD(VECTR1,VECTR2,MXI,MXF,MYI,MYF,JNDEX,XORG,YORG,       
     &   ZORG,X(1),X(2),X(3))                                           
         X(1)=X(1)-XYZ(1,II)                                            
         X(2)=X(2)-XYZ(2,II)                                            
         X(3)=X(3)-XYZ(3,II)                                            
C                                                                       
C CALCULATE HARMONIC PART OF PDF. IF THIS FACTOR IS ABOUT ZERO NO AN-   
C HARMONIC MODIFICATIONS ARE NECESSARY.                                 
C        HARMONIC PART IS     EXP(-0.5*Q(IK)*U(I)*U(K))                 
C                                                                       
 2103    EXPO  = Q(1,1)*X(1)**2 + Q(2,2)*X(2)**2 + Q(3,3)*X(3)**2 +     
     &  2.*(  Q(1,2)*X(1)*X(2) + Q(1,3)*X(1)*X(3) + Q(2,3)*X(2)*X(3) )  
C                                                                       
         IF(EXPO.LT.84.)GOTO 2553                                       
C                                                                       
C---- TO LARGE NEGATIVE EXPONENT OF EXP-FUNCTION. SET RESULT TO ZERO.   
C                                                                       
         RO=0.                                                          
         GOTO 5000                                                      
C                                                                       
C---- EXPONENTIAL FACTOR IS NOT ZERO.                                   
C                                                                       
 2553    EXPO = -EXPO*0.5                                               
         RO=EXP(EXPO)                                                   
C                                                                       
C---- ANHARMONICITY TO BE CALCULATED?                                   
C                                                                       
         IF(ITOT.NE.0) GOTO 2004                                        
C                                                                       
C---- NO ANHARMONIC PARTS DESIRED.  FINISH ONLY THE HARMONIC RESULT.    
C                                                                       
         RO=QDET*RO                                                     
         GOTO 5000                                                      
C                                                                       
C---- ANHARMONIC PART NECESSARY. CLEAR VARIABLES OF CUMULANT PRODUCTS.  
C                                                                       
 2004   CUMEX = 0.0                                                     
        CUM3  = 0.0                                                     
        CUM4  = 0.0                                                     
        CUM5  = 0.0                                                     
        CUM6  = 0.0                                                     
C                                                                       
C----  CALCULATE THE PRODUCT     Q(I,K)*X(K)   (EINSTEIN-CONVENTION)    
C                                                                       
        W(1)=Q(1,1)*X(1) + Q(1,2)*X(2) + Q(1,3)*X(3)                    
        W(2)=Q(2,1)*X(1) + Q(2,2)*X(2) + Q(2,3)*X(3)                    
        W(3)=Q(3,1)*X(1) + Q(3,2)*X(2) + Q(3,3)*X(3)                    
C                                                                       
C                                                                       
C----  IN CASE OF ALPHA FORMALISM (GC) GOTO LABEL 4000. REMEMBER THAT   
C      FOR THIS STRUCTURE FACTOR MODEL NO HERMIT POLYNOMS  ARE NECES-   
C      SARY. (OF CORSE THERE WILL BE A SELECTION IN  HERMITE POLYNOMS   
C      AND A REST FUNCTION, BUT FOR COMPUTATION IT  SHOULD  BE BETTER   
C      NOT TO SEPERATE THE COMPLETE FUNCTION.                           
C                                                                       
       IF(MODELL.EQ.4) GOTO 4000                                        
C                                                                       
C---- CALCULATE HERMITE POLYNOMS OF FOURTH ORDER AND CUMULANT PROD.     
C                                                                       
C     FROM NOW ON ONLY EDGEWORTH OR GRAM-CHARLIER POSSIBLE.             
C                                                                       
                         CALL HERMIT                                    
     &  (W,Q,ITHREE,JTHREE,CR,IFOUR,JFOUR,DR,II,CUM3,CUM4,HJKLM)        
C                                                                       
C---- CUM3 IS THIRD RANK CUMULANT PRODUCT , I.E.  C(IJK)*H(IJK)         
C                                                                       
         IF(IAWAY.EQ.1) CUM3 = 0.                                       
C                                                                       
C---- EXTENDED MAP TO BE CALCULATED?                                    
C                                                                       
 2009    IF(IEXT.EQ.0) GOTO 2112                                        
C                                                                       
C---- YES|                                                              
C                                                                       
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++      
C SECOND PART OF EXTENDED MAP  SECOND PART OF EXTENDED MAP  SECOND      
C                                                                       
C---- CALCULATE   L1(I) * HERMIT1(I)                                    
C      AT WHICH   L1(I)= X(I) - X.HARMONIC(I)                           
C                                                                       
         CUM1=0.0                                                       
         IF(JUMP2.EQ.0) GOTO 2069                                       
         DO 2007 I=1,3                                                  
 2007    CUM1=CUM1+R1(I)*W(I)                                           
C                                                                       
C---- FOURTH RANK CUMULANT PRODUCT AVAILABLE?                           
C                                                                       
 2069    CUM2 = 0.0                                                     
         CUM13= 0.0                                                     
         IF(JUMP3.EQ.0)  GOTO 2111                                      
C                                                                       
C---- YES. CALCULATE 1/2*( L1(J)*L1(K) + L2(JK) ) * HERMIT2(J,K)        
C     AT WHICH       L2(J,K)=( B.TENSOR(J,K) - B.HARMONIC(J,K) )/F      
C     AND            F = 2.PI**2                                        
C                                                                       
         DO 2070 ID=1,3                                                 
         IDD=ID                                                         
C---- IDD JUST TO PREVENT COMPILER WARNINGS.                            
 2070    CUM2=0.5*(R1(IDD)**2+R2(IDD))*HJK(IDD)+CUM2                    
         CUM2=CUM2+(R1(1)*R1(2)+R2(4))*HJK(4)                           
         CUM2=CUM2+(R1(1)*R1(3)+R2(5))*HJK(5)+(R1(2)*R1(3)+R2(6))*      
     &   HJK(6)                                                         
C                                                                       
C---- JUMP1=0 MEANS CUM13 HAS TO BE ZERO.                               
C                                                                       
         IF(JUMP1.EQ.0)  GOTO 2111                                      
C                                                                       
C---- CALCULATE THE PRODUCT L1(I)*C(JKL)*H(IJKL) (EINSTEIN CONVENTION)  
C                                                                       
                                                                        
         CUM13=4*(HJKLM(1)*R1(1)*CR(1,II) + HJKLM(2)*R1(2)*CR(2,II) +   
     &            HJKLM(3)*R1(3)*CR(3,II))                              
         CUM13=CUM13 +HJKLM(4)* ( R1(1)*CR(4,II)+CR(1,II)*R1(2) )       
         CUM13=CUM13 +HJKLM(5)* ( R1(1)*CR(6,II)+CR(1,II)*R1(2) )       
         CUM13=CUM13 +HJKLM(6)* ( R1(1)*CR(5,II)+CR(4,II)*R1(2) )*ZWIDRI
         CUM13=CUM13 +HJKLM(7)* ( R1(1)*CR(10,II)+CR(6,II)*R1(2)+       
     &            CR(4,II)*R1(3) )*EINDRI                               
         CUM13=CUM13 +HJKLM(8)* ( R1(1)*CR(7,II)+CR(6,II)*R1(3) )*ZWIDRI
         CUM13=CUM13 +HJKLM(9)* ( R1(1)*CR(2,II)+CR(5,II)*R1(2) )       
         CUM13=CUM13+HJKLM(10)* ( R1(1)*CR(8,II)+CR(10,II)*R1(2)+       
     &            CR(5,II)*R1(3) )*EINDRI                               
         CUM13=CUM13+HJKLM(11)* ( R1(1)*CR(9,II)+CR(10,II)*R1(3)+       
     &            CR(7,II)*R1(2) )*EINDRI                               
         CUM13=CUM13+HJKLM(12)* ( R1(1)*CR(3,II)+CR(7,II)*R1(3) )       
         CUM13=CUM13+HJKLM(13)* ( R1(2)*CR(8,II)+CR(2,II)*R1(3) )       
         CUM13=CUM13+HJKLM(14)* ( R1(2)*CR(9,II)+CR(8,II)*R1(3) )*ZWIDRI
         CUM13=CUM13+HJKLM(15)* ( R1(2)*CR(3,II)+CR(9,II)*R1(3) )       
         CUM13=CUM13*1.E-3                                              
C                                                                       
C---- CALCULATE CONTRIBUTION OF THE EXTENSION OF EDGEWORTH EXPANSION.   
C                                                                       
 2111    CUMEX= CUM1 + CUM2 + CUM13/TWOPI3                              
C                                                                       
C END OF PART OF EXPANSION   END OF PART OF EXPANSION   END OF PART OF  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
C                                                                       
C SAME FOR NORMAL AND EXTENDED MAPS.                                    
C                                                                       
C---- CUM4ISTHE FOURTH  CUMULANT PRODUCT.   D(IJKL)*H(IJKL)             
C                                                                       
 2112    CONTINUE                                                       
C                                                                       
 2589 IF(NO5.EQ.0) GOTO 2010                                            
C                                                                       
C ---- FIFTH RANK CUMULANT PRODUCT IS NOT ZERO. CALCULATE IT. --------  
C                                                                       
         CUM5 = RANK5(ER,JFIVE,II)                                      
C                                                                       
C ---- SIXTH ORDER CUMULANT PRODUCT IS NOT ZERO? -----------------------
C                                                                       
 2010  IF(NO6.EQ.0) GOTO 999                                            
C                                                                       
C----  YES IT IS| CALCULATE IT.                                         
C                                                                       
       CUM6 = RANK6(FR,JSIX,NO5,II)                                     
C                                                                       
C---- ACCORDING TO THE STRUCTUR FACTOR MODEL COMPUTE BRANCH.            
C                                                                       
  999 GOTO(1000,2000,3000),MODELL                                       
C                                                                       
C---- EDGEWORTH-SERIES EXPANSION.-------                                
C                                                                       
 1000 CCH=0.0                                                           
C                                                                       
C---- CALCULATE THE PRODUCT  C(I,J,K)*C(Q,P,R)*H(I,J,K,Q,P,R)           
C                                                                       
      IF(NO4.EQ.0 .OR.NO3.EQ.0) GOTO 1001                               
C                                                                       
      IF(ICCH .EQ.0) CALL   PRODUK(CCH)                                 
      IF(ICCH1.EQ.1) CUM4 = 0.                                          
C                                                                       
C---- SUMMATION OF ALL CUMULANTS.                                       
C                                                                       
 1001 RO=QDET*RO*(ISKEW  + CUMEX  + CUM3/TWOPI3 + CUM4/TWOPI4 +         
     &  0.5*CCH/TWOPI6)                                                 
C                                                                       
      GOTO 5000                                                         
C                                                                       
C                                                                       
C---- GRAM-CHARLIER-SERIES EXPANSION. -----                             
C                                                                       
 2000 RO=QDET*RO*(ISKEW +CUM3/6.+CUM4/24.+CUM5/120.+CUM6/720.)          
      GOTO 5000                                                         
C                                                                       
C---- MODIF. EDGEWORTH SERIES EXPANSION. -----                          
C                                                                       
 3000 GOTO 5000                                                         
C                                                                       
C---- MODIF. GRAM-CHARLIER SERIES EXPANSION.-----                       
C     NOW BACK TO THE ALPHA FORMALISM (GC)                              
C                                                                       
 4000 SFACT = 0.0                                                       
                         CALL JOHN                                      
     &(W,Q,ITHREE,JTHREE,CR,IFOUR,JFOUR,DR,ALF(II),CUM3,CUM4,           
     &                     SFACT)                                       
      RO=QDET*RO*(ISKEW+CUM3/6.+CUM4/24.)/(1.+SFACT)                    
      GOTO 5000                                                         
C                                                                       
C---- MULTIPLY WITH EXXAGERATION FACTOR.                                
C                                                                       
 5000    RO=EXXAG*RO                                                    
C                                                                       
C---- COUNTING OF COMPUTED POINTS AND STORING OF ONE SECTION.           
C                                                                       
         IZAEHL=IZAEHL+1                                                
         ARRAY(IZAEHL)=RO
         if(imc.eq.0) qs(izaehl) = ro                                               
C                                                                       
C---- CALCULATE MINIMUM AND MAXIMUM OF THE PROBABILITY FUNCTION.        
C                                                                       
         IF(RO.GT.ROMAX)ROMAX=RO                                        
         IF(RO.LT.ROMIN)ROMIN=RO                                        
C                                                                       
 2100    CONTINUE                                                       
 2200   CONTINUE                                                        
C                                                                       
C---- SAVE MAXIMUM OF PDF OF FIRST SECTION                              
C
c
c---- mc bit
c
      if(mc.eq.1) goto 778 
      if(imc.ne.0) call mocad(izaehl,isq)
      if(imc.eq.mc-1) call mocen(izaehl,isq)
      if(jpri.eq.0) goto 2340
 2237 continue
      goto (2231,2232,2233,2234,2235,2236) jpri
 2231 if (isq.eq.1) goto 2399
      do 2331 i=1,izaehl
 2331 array(i) = qm(i)
      jdruck = jdr(2)
      goto 2340
 2232 if (isq.eq.1) goto 2399
      do 2332 i=1,izaehl
 2332 array(i) = iminus(i)
      goto 2340
 2233 do 2333 i=1,izaehl
 2333 array(i) = qp(i)
      jdruck = jdr(3)
      goto 2340
 2234 do 2334 i=1,izaehl
 2334 array(i) = iplus(i)
      goto 2340
 2235 do 2335 i=1,izaehl
      if(qs(i).lt.0.) goto 2435
      if (isq.eq.1) qm(i) = -qp(i)
      array(i) = qs(i) + qm(i)
      if(array(i).lt.0.) array(i) = 0.
      goto 2335
 2435 array(i) = qs(i) + qp(i)
      if(array(i).gt.0.) array(i) = 0.
 2335 continue
      jdruck = jdr(4)
      goto 2340
 2236 do 2336 i=1,izaehl
      if(qs(i).lt.0.) goto 2436
      array(i) = qs(i) + sifa*qm(i)
      if(array(i).lt.0.) array(i) = 0.
      goto 2336
 2436 array(i) = qs(i) + sifa*qp(i)
      if(array(i).gt.0.) array(i) = 0.
 2336 continue
      jdruck = jdr(5)
 2340 romin = 100000.
      romax = -100000.
      do 2540 i=1,izaehl
      if(array(i).gt.romax) romax = array(i)
      if(array(i).lt.romin) romin = array(i)
 2540 continue
      if(imc.ne.0.and.imc.ne.mc-1) goto 2400                                                                   
      if(imc.eq.mc-1.and.jpri.eq.0) goto 2399
  778 IF(JZSEC-1) 777,777,781                                           
  777 IF(ROMAX.GT.0.) GOTO 780                                          
      IPOTET = 0                                                        
      INORM  = 1                                                        
      GOTO 781                                                          
  780 ROMAXT = ROMAX                                                    
      IF(INORM) 781,781,779                                             
  779 ROMAXT = ABS(PNORM)                                               
C---- PRINT SECTION AS A POINT-ARRAY, IF NPNT.NE.0 .                    
C                                                                       
  781 IF(NPNT.NE.0) GOTO 881                                            
      if(jpri.eq.2.or.jpri.eq.4) goto 881
C                                                                       
C---- DIVIDE SECTION IN STRIPS OF A MAXIMUM OF 30 COLUMNS.              
C                                                                       
      DO 880 ISPALT=1,IXA,30                                            
      IXX=ISPALT+IXI-1                                                  
      NXA=MIN0(30,IXF-IXX+1)                                            
      NX=0                                                              
C                                                                       
C---- PRINT HEADING OF SECTION.                                         
C                                                                       
            DO 782 I=1,NXA,2                                            
            NX=NX+1                                                     
  782       FX(NX)=I+IXX-1                                              
      J=JZSEC+IZI-1
      if(imc.eq.0)jsave=j                                            
      if(imc.ne.0)j=jsave 
      IF(IPOTET.EQ.0)                                                   
     &WRITE(IOUT,8)(TITLE(L),L=1,17),jdruck,J,IZ,(TAG(J,II),J=1,2),ROMIN
     &,ROMAX,(FX(I),I=1,NX)                                                   
    8 FORMAT(1H1,17A4,17X,a25,4x,7HSECTION,I4,1H/,I4,/,
     &/                                        38X,41H PROBABILITY-MAP O
     &F THE ATOM LABELED       ,2A4,//32X,41H FLUCTUATION OF THE PDF IN 
     &THIS SECTION       ,F8.0,4X,3H-  ,F8.0//F10.0,14F8.0,//)          
      IF(IPOTET.EQ.1)                                                   
     &WRITE(IOUT,9)(TITLE(L),L=1,17),J,IZ,(TAG(J,II),J=1,2),T,ROMIN,    
     & ROMAX,(FX(I),I=1,NX)                                             
    9 FORMAT(1H1,17A4,29X,7HSECTION,I4,1H/,I4//38X,37H POTENTIAL - MAP O
     &F THE ATOM LABELED ,4X,2A4,//38X,28H MEASURING  TEMPERATURE  WAS, 
     &F8.0,2X,6HKELVIN //38X,35H OUTPUT SHOWS THE POTENTIAL IN MEV.     
     &///32X,43H FLUCTUATION OF THE PDF IN THIS SECTION WAS             
     & ,F8.0,4X,3H-  ,F8.0//F10.0,14F8.0,//)                            
C                                                                       
                                                                        
      DO 880 JYSEC=1,IYA                                                
      IVORN=(JYSEC-1)*IXA                                               
      FY=JYSEC+IYI-1                                                    
      L=0                                                               
            DO 870 JXSEC=ISPALT,ISPALT+NXA-1                            
            IPSV=IVORN+JXSEC-1                                          
            JXX =JXSEC-ISPALT+1                                         
  805       IPSV=IPSV+1                                                 
C                                                                       
C----- CALCULATE POTENTIAL OUT OF PDF. RESTRICT IT TO SENSFUL VALUES.   
C      IF PDF IS NEGATIVE MARK THE POTENTIAL WITH A SIGN.               
C                                                                       
               IF (IPOTET) 806,806,807                                  
 807           RO=ARRAY(IPSV)                                           
               RO=RO/ROMAXT                                             
               IF(RO.LT.1.E-20) GOTO 808                                
               ARRAY(IPSV) = -0.08617*T*ALOG(ABS(RO))                   
               GOTO 806                                                 
C----          INSTEAD OF EXACTLY ZERO 1.E-25 SHOULD BEEN WRITTEN       
C              TO PREVENT ERRORS OF THE PLOTTING PROGRAMME.             
 808           ARRAY(IPSV)=1.E-25                                       
C                                                                       
 806        RO=ARRAY(IPSV)                                              
C                                                                       
C---- CALCULATE OUTPUT AND OUTPUT FORMAT FOR EACH LINE.                 
C                                                                       
            J=RO+SIGN(.5,RO)                                            
            IF(J) 810,860,830                                           
  810       IF(J+1000)820,820,850                                       
  820       J=J+1000                                                    
            GO TO 810                                                   
  830       IF(J-10000)850,840,840                                      
  840       J=J-10000                                                   
            GO TO 830                                                   
  850       IF(J.GT.IMINI.AND.J.LT.IMAXI) GOTO 860                      
            L=L+1                                                       
            IP(L)=J                                                     
            FRMT(JXX + 3)=NZRO                                          
            GO TO 870                                                   
  860       FRMT(JXX + 3)=ZRO                                           
  870       CONTINUE                                                    
  876 FRMT(NXA+4)=ENDD                                                  
C                                                                       
C---- IN CORRESPONDENCE TO OUTPUT ANGLE BRANCHING.                      
C                                                                       
      IF(IS)1021,1022,1023                                              
 1021 LSTORE=L                                                          
      ICHECK=(JYSEC-1)*IHILF+1                                          
      FRMT(3)=IU(ICHECK)                                                
      CHEKK = FLOAT(ICHECK/4)
      ICHECK=INT(CHEKK)+1                                            
      IF(ICHECK+NXA-30)1043,1043,1042                                   
 1042 DO 1071 IO=34-ICHECK,3+NXA                                        
      IF(FRMT(IO).EQ.NZRO)L=L-1                                         
 1071 FRMT(IO)=LEER                                                     
 1043 IF(L.EQ.0)GOTO 872                                                
      WRITE(IOUT,FRMT)LSPC,FY,(IP(I),I=1,L)                             
      GO TO 880                                                         
 1022 FRMT(3)=IU(1)                                                     
      GOTO 877                                                          
 1023 IF(JYSEC.EQ.1)GO TO 1050                                          
      IF(IHILF.EQ.3) GOTO 1201                                          
      IF(LUECK.EQ.0)ISTART=ISTART+1                                     
      LUECK=MOD(4+LUECK-IHILF,4)                                        
      FRMT(3)=IU(1+LUECK)                                               
      ICOR=0                                                            
      IF(ISTART.GT.NXA)GOTO 880                                         
      DO 1054 IKLJ=1,ISTART-1                                           
      IF(FRMT(3+IKLJ).EQ.NZRO)ICOR=ICOR+1                               
 1054 FRMT(3+IKLJ)=LEER                                                 
      ISTAR=ICOR+1                                                      
 1045 IF(L.EQ.0)GOTO 872                                                
      WRITE(IOUT,FRMT)LSPC,FY,(IP(I),I=ISTAR,L)                         
      GO TO 880                                                         
 1201 ISTART=ISTART+1                                                   
      IF(LUECK.EQ.0) LUECK=4                                            
      LUECK=LUECK-1                                                     
      IF(ISTART.EQ.NXA) GOTO 872                                        
      GOTO 1045                                                         
 1050 ISTART=1                                                          
      LUECK=0                                                           
      ISTAR=1                                                           
      FRMT(3)=IU(1)                                                     
      GOTO 1045                                                         
 877  IF(L.EQ.0) GO TO 872                                              
      WRITE(IOUT,FRMT)LSPC,FY,(IP(I),I=1,L)
      GOTO 880
C
C---- Output of a line which elemnts are all zero
C
 872  WRITE(IOUT,6) LSPC,FY
   6  FORMAT(A1,F4.0)                             
 880  CONTINUE                                                       
C                                                                       
C----  SECTION NOW COMPLETE - DUMP PLOT  OUTPUT IF NEEDED               
C                                                                       
 881  IF(IPLT.EQ.0)GO TO 890                                            
      J=JZSEC+IZI-1                                                     
C---- CHANGE OUTPUT FOR SPECIAL PLANE                                   
C                                                                       
      IF(IPLANE.EQ.0) GOTO 809                                          
      J=1                                                               
  809 GOTO (8091,8092,890,8093,890,8094,8095) JPRI + 1
 8091 IPLOTT = S2
      GOTO 8090
 8092 IF(ISQ.EQ.1)GOTO 890
      IPLOTT = S3
      GOTO 8090
 8093 IPLOTT = S4
      GOTO 8090
 8094 IPLOTT = S5
      GOTO 8090
 8095 IPLOTT = S6    
 8090 WRITE(IPLOTT)IXA,IYA,J,IZ,IX,IXI,IXF,IY,IYI,IYF,IZI,IZF,IATOM         
      WRITE(IPLOTT)(ARRAY(I),I=1,IZAEHL)
c     type 8089, iplott
c8089 format(' FILE FOR0',i2,' WRITTEN.')                                                                                              
C                                                                       
C---- RESTORE EXTENSION PARAMETER                                       
C                                                                       
      IEXT=IHOLD
      if(mc.eq.1) goto 2400
 2399 if(imc.eq.mc-1) jpri = jpri + 1
      if(imc.eq.mc-1.and.jpri.lt.7) goto 2237 
  890 CONTINUE                         
c
c---- add monte carlo sections
c
 2400 CONTINUE
C                                                                       
 2500 continue                                                          
C                                                                       
C        END OF MAP.                                                    
C                                                                       
  900 RETURN                                                            
C                                                                       
C---- SOME AUXILLARY PARTS.                                             
C                                                                                                                                
C                                                                       
C---- ERRORMESSAGES, WHICH CAUSES CHANGES IN PROGRAMME FLOW.            
C                                                                       
  908 WRITE(IOUT,909) MODELL                                            
  909 FORMAT(////13H0ERRORMESSAGE/1X,12(1H-)//23H FOR THE GIVEN STRUCTUR
     &,20HE FACTOR MODEL  NO.,I4,2X,30H IS NO EXTENDED MAP AVAILABLE.// 
     &27H0INPUT RESET TO NORMAL MAP.//1H )                              
      GOTO 241                                                          
C                                                                       
  910 WRITE(IOUT,911)                                                   
  911 FORMAT( 4(/),13H0ERRORMESSAGE/1X,12H(1H-)//                       
     &        85H0YOU HAVE TO INCREASE THE DIMENSION OF THE ARRAY *ARRAY
     &* TO THE ABOVE PRINTED VALUE.   //32H0I'M SORRY, I'VE TO CALL  EXI
     &   ,2X,1H.)                                                       
      CALL EXIT                                                         
C                                                                       
  912 WRITE(IOUT,913)MODELL                                             
  913 FORMAT(4(/),13H0ERRORMESSAGE/1X,12(1H-)//19H FOR THE STRUCTURE ,  
     & 5HMODEL,I4,2X,38HTHERE IS NO PDF-CALCULATION AVAILABLE.//        
     & 31HI'M SORRY, I'VE TO CALL 'EXIT'.)                              
      CALL EXIT                                                         
C                                                                       
  917 WRITE(IOUT,918)MODELL                                             
  918 FORMAT(4(/),13H0ERRORMESSAGE/1X,12(1H-)//                         
     & 31H FOR THE STRUCTURE-FACTOR MODEL ,I4,2X,18HTHERE IS NO FIFTH , 
     & 26HOR HIGHER ORDER AVAILABLE.//29H SELECTED TYPE WILL BE RESET.  
     & //1H )                                                           
      IPDFF(4)=0                                                        
      IPDFF(5)=0                                                        
      GOTO 242                                                          
C                                                                       
  612 WRITE(IOUT,919)                                                   
  919 FORMAT(4(/),13HERRORMESSAGE/1X,12(1H-)//21H NOT ANY CUMULANT HAS  
     & 15H BEEN SELECTED.//22H RESET TO 'TOTAL' MAP.//1H )              
      DO 4139 I=1,5                                                     
 4139 IPDFF(I)=1                                                        
      GOTO 414                                                          
C                                                                       
  921 WRITE(IOUT,922)                                                   
  922 FORMAT(4(/),13H0ERRORMESSAGE/1X,12(1H-)//19H THE SELECTED TENSO,  
     &  48HR(S) IS (ARE) NULL. SO THIS MESSAGE IS A RESULT.//           
     & 28H THAT'S ALL. JOB TERMINATED.)                                 
      CALL EXIT                                                         
C                                                                       
  923 FORMAT(4(/),13H0ERRORMESSAGE/1X,12(1H-)//19H CALCULATION OF POT,  
     &  54HENTIAL WITHOUT  SECOND ORDER CUMULANTS IS MEANINGLESS./      
     &  54H PROGRAMME ADD THE SECOND ORDER CUMULANT COEFFICIENTS.//1H ) 
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE ACHSE(ISECT,IX,IXI,IXF,IY,IYI,IYF,IZ,IZI,IZF)          
C                                                                       
C     TRANSFORMATION OF GRIDS, WHICH IS ALWAYS  X,Y,Z DURING INPUT      
C     TO THAT SEQUENCE, WHICH IS NECESSARY FOR SELECTED SECTIONS.       
C                                                                       
C                            BY UDO ZUCKER                              
C                                                                       
      IF(ISECT)25,25,5                                                  
    5 GO TO (10,30,50,60,70),ISECT                                      
C                                                                       
C---- XZY SECTIONS                                                      
C                                                                       
   10 IW  = IY                                                          
      IWI = IYI                                                         
      IWF = IYF                                                         
   20 IY =IZ                                                            
      IYI=IZI                                                           
      IYF=IZF                                                           
      IZ =IW                                                            
      IZI=IWI                                                           
      IZF=IWF                                                           
   25 RETURN                                                            
C                                                                       
C---- YXZ SECTIONS                                                      
C                                                                       
   30 IW = IY                                                           
      IWI=IYI                                                           
      IWF=IYF                                                           
   40 IY=IX                                                             
      IYI=IXI                                                           
      IYF=IXF                                                           
      IX=IW                                                             
      IXI=IWI                                                           
      IXF=IWF                                                           
      RETURN                                                            
C                                                                       
C---- YZX SECTIONS                                                      
C                                                                       
   50 IW=IX                                                             
      IWI=IXI                                                           
      IWF=IXF                                                           
      IX=IY                                                             
      IXI=IYI                                                           
      IXF=IYF                                                           
      GO TO 20                                                          
C                                                                       
C---- ZYX SECTIONS                                                      
C                                                                       
   60 IW=IX                                                             
      IWI=IXI                                                           
      IWF=IXF                                                           
      IX=IZ                                                             
      IXI=IZI                                                           
      IXF=IZF                                                           
      IZ=IW                                                             
      IZI=IWI                                                           
      IZF=IWF                                                           
      RETURN                                                            
C                                                                       
C---- ZXY SECTIONS                                                      
C                                                                       
   70 IW=IZ                                                             
      IWI=IZI                                                           
      IWF=IZF                                                           
      IZ=IY                                                             
      IZI=IYI                                                           
      IZF=IYF                                                           
      GO TO 40                                                          
C                                                                       
      END                                                               
C                                                                       
                          SUBROUTINE HERMIT                             
     &   (W,Q,ITHREE,JTHREE,CR,IFOUR,JFOUR,DR,II,CUM3,CUM4,HJKLM)       
C                                                                       
C     THIS ROUTINE CALCULATES THE  PRODUCTS OF  HERMITE  POLYNOMS  OF   
C     THIRD AND FOURTH ORDER BY  CUMULANTS  OF  CORRESPONDING  ORDER.   
C     NATURATELLY ONE COULD USE THE  SUBROUTINE  JOHN WITH ALF=0, BUT   
C     THIS ROUTINE HERE IS FASTER.                                      
C                                                                       
C                                                                       
C                    UDO ZUCKER  NOVEMBER 1980                          
C                                                                       
C     IN THE ROUTINE 3 QUANTITYS WILL BE  CALCULATED  CUM3, CUM4 AND    
C     THE ARRAY HJKLM HOLDING THE HERMITE POLYNOMS  OF FOURTH ORDER.    
C                                                                       
      DIMENSION W(3),Q(3,3),CR(10,60),DR(15,60),JTHREE(10),JFOUR(15)    
     & ,HJKLM(15)                                                       
C                                                                       
      W1    =  W(1)                                                     
      W2    =  W(2)                                                     
      W3    =  W(3)                                                     
      W11   = W1*W1                                                     
      W22   = W2*W2                                                     
      W33   = W3*W3                                                     
C                                                                       
C                                                                       
C---- NO THIRD ORDER PRODUCT,IF ITHREE EQUAL ZERO.                      
C                                                                       
      IF(ITHREE.EQ.0) GOTO 3000                                         
C                                                                       
C---- INITIALIZE                                                        
C                                                                       
      TERM  = 0.                                                        
      TERMZ = 0.                                                        
C                                                                       
      IF(JTHREE(1))1,2,1                                                
    1 TERM1 = W11*W1                                                    
      TERM2 = 3.*W1*Q(1,1)                                              
      TERM  = (TERM1-TERM2)*CR(1,II)                                    
    2 IF(JTHREE(2))3,4,3                                                
    3 TERM1 = W22*W2                                                    
      TERM2 = 3.*W2*Q(2,2)                                              
      TERM  = TERM + (TERM1-TERM2)*CR(2,II)                             
    4 IF(JTHREE(3))5,6,5                                                
    5 TERM1 = W33*W3                                                    
      TERM2 = 3.*W3*Q(3,3)                                              
      TERM  = TERM + (TERM1-TERM2)*CR(3,II)                             
    6 IF(JTHREE(4))7,8,7                                                
    7 TERM1 = W11*W2                                                    
      TERM2 = (2.*W1*Q(1,2)+W2*Q(1,1))                                  
      TERMZ = (TERM1-TERM2)*CR(4,II)                                    
    8 IF(JTHREE(5))9,10,9                                               
    9 TERM1 = W22*W1                                                    
      TERM2 = (W1*Q(2,2)+2.*W2*Q(1,2))                                  
      TERMZ = TERMZ + (TERM1-TERM2)*CR(5,II)                            
   10 IF(JTHREE(6))11,12,11                                             
   11 TERM1 = W11*W3                                                    
      TERM2 = (2.*W1*Q(1,3)+W3*Q(1,1))                                  
      TERMZ = TERMZ + (TERM1-TERM2)*CR(6,II)                            
   12 IF(JTHREE(7))13,14,13                                             
   13 TERM1 = W1*W33                                                    
      TERM2 = (W1*Q(3,3)+2.*W3*Q(1,3))                                  
      TERMZ = TERMZ + (TERM1-TERM2)*CR(7,II)                            
   14 IF(JTHREE(8))15,16,15                                             
   15 TERM1 = W22*W3                                                    
      TERM2 = (2.*W2*Q(2,3)+W3*Q(2,2))                                  
      TERMZ = TERMZ + (TERM1-TERM2)*CR(8,II)                            
   16 IF(JTHREE(9))17,18,17                                             
   17 TERM1 = W2*W33                                                    
      TERM2 = (W2*Q(3,3)+2.*W3*Q(2,3))                                  
      TERMZ = TERMZ + (TERM1-TERM2)*CR(9,II)                            
   18 TERM  = TERM + TERMZ*3.0                                          
      TERMZ = 0.0                                                       
      IF(JTHREE(10))19,20,19                                            
   19 TERM1 = W1*W2*W3                                                  
      TERM2 = (W1*Q(2,3)+W2*Q(1,3)+W3*Q(1,2))                           
      TERMZ = (TERM1-TERM2)*CR(10,II)*6.0                               
   20 CUM3  = TERM + TERMZ                                              
      CUM3  = CUM3 * 1.E-3                                              
C                                                                       
C---- IF IFOUR=0 FOURTH RANK PRODUCTS ARE NOT REQUIRED.                 
C                                                                       
 3000 IF(IFOUR.EQ.0) GOTO 4000                                          
C                                                                       
      TERM  = 0.                                                        
      TERMZ = 0.                                                        
C                                                                       
      IF(JFOUR(1))101,102,101                                           
  101 TERM1    = W11*W11                                                
      TERM2    = 6.*W11*Q(1,1)                                          
      TERM3    = 3.*Q(1,1)**2                                           
      HJKLM(1) = TERM1-TERM2+TERM3                                      
      TERM     = HJKLM(1)*DR(1,II)                                      
  102 IF(JFOUR(2))103,104,103                                           
  103 TERM1    = W22*W22                                                
      TERM2    = 6.*W22*Q(2,2)                                          
      TERM3    = 3.*Q(2,2)**2                                           
      HJKLM(2) = TERM1-TERM2+TERM3                                      
      TERM     = TERM + HJKLM(2)*DR(2,II)                               
  104 IF(JFOUR(3))105,106,105                                           
  105 TERM1    = W33*W33                                                
      TERM2    = 6.*W33*Q(3,3)                                          
      TERM3    = 3.*Q(3,3)**2                                           
      HJKLM(3) = TERM1-TERM2+TERM3                                      
      TERM     = TERM + HJKLM(3)*DR(3,II)                               
  106 IF(JFOUR(4))107,108,107                                           
  107 TERM1    = W11*W1*W2                                              
      TERM2    = (W11*Q(1,2)*3.+W1*W2*Q(1,1)*3.)                        
      TERM3    = 3.*Q(1,1)*Q(1,2)                                       
      HJKLM(4) = TERM1-TERM2+TERM3                                      
      TERMZ    = HJKLM(4)*DR(4,II)                                      
      HJKLM(4) = HJKLM(4)*4.0                                           
  108 IF(JFOUR(5))109,110,109                                           
  109 TERM1    = W11*W1*W3                                              
      TERM2    = (W11*Q(1,3)*3.+W1*W3*Q(1,1)*3.)                        
      TERM3    = 3.*Q(1,1)*Q(1,3)                                       
      HJKLM(5) = TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(5)*DR(5,II)                              
      HJKLM(5) =HJKLM(5)*4.0                                            
  110 IF(JFOUR(9))111,112,111                                           
  111 TERM1    = W22*W2*W1                                              
      TERM2    = (W22*Q(1,2)*3.+W2*W1*Q(2,2)*3.)                        
      TERM3    = 3.*Q(2,2)*Q(1,2)                                       
      HJKLM(9) = TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(9)*DR(9,II)                              
      HJKLM(9) = HJKLM(9)*4.0                                           
  112 IF(JFOUR(12))113,114,113                                          
  113 TERM1    = W33*W3*W1                                              
      TERM2    = (W33*Q(1,3)*3.+W1*W3*Q(3,3)*3.)                        
      TERM3    = 3.*Q(3,3)*Q(1,3)                                       
      HJKLM(12)= TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(12)*DR(12,II)                            
      HJKLM(12)= HJKLM(12)*4.0                                          
  114 IF(JFOUR(13))115,116,115                                          
  115 TERM1    = W22*W2*W3                                              
      TERM2    = (W22*Q(2,3)*3.+W2*W3*Q(2,2)*3.)                        
      TERM3    = 3*Q(2,2)*Q(2,3)                                        
      HJKLM(13)= TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(13)*DR(13,II)                            
      HJKLM(13)= HJKLM(13)*4.0                                          
  116 IF(JFOUR(15))117,118,117                                          
  117 TERM1    = W33*W3*W2                                              
      TERM2    = (W33*Q(2,3)*3.+W2*W3*Q(3,3)*3.)                        
      TERM3    = 3.*Q(3,3)*Q(2,3)                                       
      HJKLM(15)= TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(15)*DR(15,II)                            
      HJKLM(15)= HJKLM(15)*4.0                                          
  118 TERM     = TERM + TERMZ*4.0                                       
      TERMZ    = 0.0                                                    
      IF(JFOUR(6))119,120,119                                           
  119 TERM1    = W11*W22                                                
      TERM2    = (W11*Q(2,2)+4.*W1*W2*Q(1,2)+W22*Q(1,1))                
      TERM3    = Q(1,1)*Q(2,2)+2.*Q(1,2)**2                             
      HJKLM(6) = TERM1-TERM2+TERM3                                      
      TERMZ    = HJKLM(6)*DR(6,II)                                      
  120 IF(JFOUR(8))121,122,121                                           
  121 TERM1    = W11*W33                                                
      TERM2    = (W11*Q(3,3)+4.*W1*W3*Q(1,3)+W33*Q(1,1))                
      TERM3    = Q(1,1)*Q(3,3)+2.*Q(1,3)**2                             
      HJKLM(8) = TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(8)*DR(8,II)                              
      HJKLM(8) = HJKLM(8)*6.0                                           
  122 IF(JFOUR(14))123,124,123                                          
  123 TERM1    = W33*W22                                                
      TERM2    = (W33*Q(2,2)+4.*W3*W2*Q(2,3)+W22*Q(3,3))                
      TERM3    = Q(3,3)*Q(2,2)+2.*Q(2,3)**2                             
      HJKLM(14)= TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(14)*DR(14,II)                            
      HJKLM(14)= HJKLM(14)*6.0                                          
  124 TERM     = TERM + TERMZ*6.0                                       
      TERMZ    = 0.0                                                    
      IF(JFOUR(10))125,126,125                                          
  125 TERM1    = W22*W3*W1                                              
      TERM2    = (W22*Q(1,3)+2.*W2*W3*Q(1,2)+2.*W2*W1*Q(2,3) +          
     &           W3*W1*Q(2,2))                                          
      TERM3    = Q(2,2)*Q(1,3)+2.*Q(2,3)*Q(1,2)                         
      HJKLM(10)= TERM1-TERM2+TERM3                                      
      TERMZ    = HJKLM(10)*DR(10,II)                                    
      HJKLM(10)= HJKLM(10)*12.0                                         
  126 IF(JFOUR(11))127,128,127                                          
  127 TERM1    = W1*W2*W33                                              
      TERM2    = (W1*W2*Q(3,3)+W2*W3*Q(1,3)*2.+W33*Q(1,2) +             
     &           2.*W1*W3*Q(2,3))                                       
      TERM3    = Q(1,2)*Q(3,3)+2.*Q(1,3)*Q(2,3)                         
      HJKLM(11)= TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(11)*DR(11,II)                            
      HJKLM(11)= HJKLM(11)*12.0                                         
  128 IF(JFOUR(7))129,130,129                                           
  129 TERM1    = W11*W2*W3                                              
      TERM2    = (W11*Q(2,3)+2.*W1*W2*Q(1,3)+2.*W1*W3*Q(1,2) +          
     &           W2*W3*Q(1,1))                                          
      TERM3    = Q(1,1)*Q(2,3)+2.*Q(1,2)*Q(1,3)                         
      HJKLM(7) = TERM1-TERM2+TERM3                                      
      TERMZ    = TERMZ + HJKLM(7)*DR(7,II)                              
      HJKLM(7) = HJKLM(7)*12.0                                          
  130 CUM4     = TERM + TERMZ*12.0                                      
      CUM4     = CUM4 * 1.E-4                                           
C                                                                       
C                                                                       
 4000 RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE UMKEHR(BETA,I,QDET)                                    
C                                                                       
C     INVERSION OF THE SYMETRIC TENSOR BETA.                            
C                                                                       
C                                                                       
      DIMENSION BETA(6,60)                                              
C                                                                       
      COMMON /INVBET/Q(3,3)                                             
      DATA TWOPIQ/19.739208/
C                                                                       
      Q(1,1)=BETA(2,I)*BETA(3,I)-BETA(6,I)**2                           
      Q(1,2)=-BETA(4,I)*BETA(3,I)+BETA(5,I)*BETA(6,I)                   
      Q(1,3)=BETA(4,I)*BETA(6,I)-BETA(5,I)*BETA(2,I)                    
      QDET  =BETA(1,I)*Q(1,1)+BETA(4,I)*Q(1,2)+BETA(5,I)*Q(1,3)         
      Q(1,1)=Q(1,1)/QDET                                                
      Q(1,2)=Q(1,2)/QDET                                                
      Q(1,3)=Q(1,3)/QDET                                                
      Q(2,2)=(BETA(1,I)*BETA(3,I)-BETA(5,I)**2)/QDET                    
      Q(2,3)=(-BETA(1,I)*BETA(6,I)+BETA(5,I)*BETA(4,I))/QDET            
      Q(3,3)=(BETA(1,I)*BETA(2,I)-BETA(4,I)**2)/QDET                    
      Q(2,1)=Q(1,2)                                                     
      Q(3,1)=Q(1,3)                                                     
      Q(3,2)=Q(2,3)                                                     
      DO 1 J=1,3                                                        
C                                                                       
C---- NORMALIZE TO A VALUE NEEDED IN THE PDF PROGRAM.                   
C                                                                       
         DO 1 K=1,3                                                     
   1     Q(J,K)=TWOPIQ*Q(J,K)                                           
C                                                                       
      QDET=5.568328/SQRT(QDET)
C     qdet = 1./(15.74961*sqrt(qdet)) !5.568328 is pi**1.5 we want 2pi**1.5
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PRODUK(CCH)                                            
C                                                                       
      COMMON /KUMU/C111,C222,C333,C112,C122,C113,C133,C223,C233,C123    
C                                                                       
C     THIS ROUTINE CALCULATES THE PRODUCT                               
C                                                                       
C            JKL    QRS                                                 
C     CCH = C    * C    * H          (EINSTEIN CONVENTION)              
C                          JKLQRS                                       
C                                                                       
C---- CALCULATE OFTEN NEEDED FORMULAAS.                                 
C                                                                       
      CALL START                                                        
C                                                                       
C---- CALCULATE NOW THE PRODUCT. WHAT FUN |                             
C                                                                       
      T1=C111*C111*H1(1)+C222*C222*H1(2)+C333*C333*H1(3)                
      T2=C111*(C112*H2(1,2)+C113*H2(1,3))+C222*(C122*H2(2,1)+C223*      
     &   H2(2,3))+C333*(C133*H2(3,1)+C233*H2(3,2))                      
      T3=(2.*C111*C122+3.*C112*C112)*H3(1,2)+(2.*C111*C133+3.*C113*     
     &  C113)*H3(1,3)+(2.*C222*C112+3.*C122*C122)*H3(2,1)+(2.*C222      
     &  *C233+3.*C223*C223)*H3(2,3)+(2.*C333*C113+3.*C133*C133)*H3      
     &  (3,1)+(2.*C333*C223+3.*C233*C233)*H3(3,2)                       
      T4=(3.*C113*C112+2.*C111*C123)*H4(1,2,3)+(3.*C223*C122+2.*        
     &  C222*C123)*H4(2,1,3)+(3.*C133*C233+2.*C123*C333)*H4(3,1,2)      
      T5=(C111*C222+9.*C122*C112)*H7(1,2)+(C111*C333+9.*C133*C113)      
     & *H7(1,3)+(C222*C333+9.*C223*C233)*H7(2,3)                        
      T6=(C111*C223+6.*C112*C123+3.*C122*C113)*H5(1,2,3)+(C111*C233     
     & +6.*C113*C123+3.*C112*C133)*H5(1,3,2)+(C222*C113+6.*C122*C123    
     & +3.*C112*C223)*H5(2,1,3)+(C133*C222+6.*C223*C123+3.*C122*C233    
     & )*H5(2,3,1)+(C333*C112+6.*C133*C123+3.*C113*C233)*H5(3,1,2)+     
     & (C333*C122+6.*C233*C123+3.*C223*C133)*H5(3,2,1)                  
      T7=(C112*C233+C113*C223+C133*C122+2.*C123*C123)*H6(1,2,3)         
C                                                                       
      CCH = T1 + 6.*T2 + 3.*T3 + 6.*T4 + 2.*T5 +  6.*T6 + 18.*T7        
      CCH = CCH*1.E-6
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE START                                                  
C                                                                       
      COMMON /INVBET/ Q(3,3)                                            
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
C                                                                       
      DO 1 IT=1,3                                                       
      W2(IT)=W(IT)**2                                                   
      W3(IT)=W2(IT)*W(IT)                                               
      W4(IT)=W3(IT)*W(IT)                                               
      W5(IT)=W4(IT)*W(IT)                                               
    1 CONTINUE                                                          
C                                                                       
C---- THIS TERMS ARE OFTEN USED IN CALCULATION OF HERMIT POLYN.         
C     OF HIGHER THAN FIFTH ORDER.                                       
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE PREPAR(MODELL,II,CR)                                   
C                                                                       
C     CALCULATION OF SOME OFTEN USED VALUES.                            
C                                                                       
      COMMON /INVBET/Q(3,3)                                             
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /KUMU/C111,C222,C333,C112,C122,C113,C133,C223,C233,C123    
C                                                                       
      DIMENSION CR(10,60)                                               
C                                                                       
      DO 1 I=1,3                                                        
      P(I) = Q(I,I)                                                     
      P2(I) = P(I)**2                                                   
   1  P3(I) = P2(I)*P(I)                                                
C                                                                       
      IF(MODELL.NE.1) RETURN                                            
C                                                                       
      C111  = CR(1,II)                                                  
      C222  = CR(2,II)                                                  
      C333  = CR(3,II)                                                  
      C112  = CR(4,II)                                                  
      C122  = CR(5,II)                                                  
      C113  = CR(6,II)                                                  
      C133  = CR(7,II)                                                  
      C223  = CR(8,II)                                                  
      C233  = CR(9,II)                                                  
      C123  = CR(10,II)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HH1(J)                                                   
C                                                                       
C     HERMIT POLYNOM OF FIFTH RANK   HH1=H(J,J,J,J,J)                   
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      HH1 = W(J) * ( W4(J)-10.*P(J)*( W2(J)-1.5*P(J) )  )               
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HH2(J,L)                                                 
C                                                                       
C     HERMIT POLYNOM OF FIFTH RANK   HH2=H(J,J,J,J,L)                   
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
C                                                                       
C     HH2 = W(J)*( W(L)*(W3(J)-6.*P(J)*W(J)) - 4.*Q(J,L)*               
C    & (W2(J)-3.*P(J)) )  + 3.*P2(J)*W(L)                               
      HZW = W(L)*(W3(J)-6.*P(J)*W(J))-4.*Q(J,L)*(W2(J)-3.*P(J))
      HH2 = W(J)*HZW + 3.*P2(J)*W(L)
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HH3(J,L)                                                 
C                                                                       
C     HH3 = H(J,J,J,L,L)                                                
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
C                                                                       
      HH3 = W(J)* ( W(L)*( W(L)*W2(J)-6.*Q(J,L)*W(J)-3.*P(J)*W(L) )     
     & - P(L)* (W2(J)-3.*P(J)) ) + 6.*Q(J,L)*(Q(J,L)*W(J)+P(J)*W(L))    
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HH4(I,J,L)                                               
C                                                                       
C     HERMIT POLYNOM OF FIFTH RANK   HH4=H(I,I,I,J,L)                   
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      HH4 = W(I)*( W2(I)*W(J)*W(L)-W(I)*(3.*Q(I,L)*W(J)+3.*Q(I,J)*      
     & W(L)+Q(J,L)*W(I))  - 3.*(P(I)*(W(J)*W(L)-Q(J,L))-2.*Q(I,J)*      
     & Q(I,L))) + 3.*P(I)*(Q(I,L)*W(J)+Q(I,J)*W(L))                     
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HH5(I,J,L)                                               
C                                                                       
C     HERMIT POLYNOM OF FIFTH RANK   HH5 = H(I,I,J,J,L)                 
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      HH5 = W(I)* (W(J)* (W(J)* (W(I)*W(L)-2.*Q(I,L)) -2.*Q(J,L)*W(I)   
     & -4.*Q(I,J)*W(L)) -P(J)*W(I)*W(L) + 2.*P(J)*Q(I,L)+4.*Q(I,J)*     
     & Q(J,L))  + W(J)* (4.*Q(I,J)*Q(I,L)-P(I)*(W(J)*W(L)-2.*Q(J,L)))   
     & + W(L)* (2.*Q(I,J)**2 + P(I)*P(J) )                              
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H1(I)                                                    
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK     H1=H(I,I,I,I,I,I)                
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      H1=W4(I)* (W2(I)-15.*P(I))  +45.*W2(I)*P2(I)-15.*P3(I)            
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H2(I,J)                                                  
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK     H2=H(I,I,I,I,I,J)                
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
C     H2=W3(I)*(W2(I)*W(J)-5.*W(I)*Q(I,J)-10.*W(J)*P(I))+15.*W(I)       
C    & *W(J)*P2(I)+30.*W2(I)*Q(I,J)*P(I)-15.*Q(I,J)*P2(I)               
      ZW = W2(I)*W(J)-5.*W(I)*Q(I,J)-10.*W(J)*P(I)
      H2 = W3(I)*ZW + 15.*W(I)*W(J)*P2(I)+30.*W2(I)*Q(I,J)*P(I)
     &- 15.*Q(I,J)*P2(I)
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H3(I,J)                                                  
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK     H3 = H(I,I,I,I,J,J)              
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
                                                                        
      H3=W4(I)*W2(J)-W4(I)*P(J)-8.*W3(I)*W(J)*Q(I,J)-6.*W2(I)*W2(J)*    
     & P(I)+6.*W2(I)*P(I)*P(J)+12.*W2(I)*Q(I,J)*Q(I,J)+24.*W(I)*W(J)    
     & *Q(I,J)*P(I)+3.*W2(J)*P2(I)-3.*P2(I)*P(J)-12.*Q(I,J)*Q(I,J)*P    
     & (I)                                                              
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H4(I,J,L)                                                
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK     H4=H(I,I,I,I,J,L)                
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      H4= W4(I)*W(J)*W(L)-W4(I)*Q(J,L)-4.*W3(I)*W(L)*Q(I,J)-4.*W3(I)*   
     & W(J)*Q(I,L)-6.*W2(I)*W(J)*W(L)*P(I)+6.*W2(I)*P(I)*Q(J,L)+12.*    
     & W2(I)*Q(I,L)*Q(I,J)+12.*W(I)*W(J)*P(I)*Q(I,L)+12.*W(I)*W(L)*P(   
     & I)*Q(I,J)+3.*W(J)*W(L)*P(I)*P(I)-12.*P(I)*P(I)*Q(J,L)-3.*P(I)*   
     & Q(I,L)*P(J)                                                      
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H5(I,J,L)                                                
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK      H5 = H(I,I,I,J,J,L)             
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      T1= W3(I)*W2(J)*W(L)-2.*W3(I)*W(J)*Q(J,L)-W3(I)*W(L)*P(J)-3.*W2(  
     &I)*W2(J)*Q(I,L)-6.*W2(I)*W(J)*W(L)*Q(I,J)-3.*W(I)*W2(J)*W(L)*P(I  
     & )+6.*W2(I)*Q(J,L)*Q(I,J)+3.*W2(I)*Q(I,L)*P(J)+6.*W(I)*W(J)*P(I)  
     & *Q(J,L)+12.*W(I)*W(J)*Q(I,J)*Q(I,L)+6.*W(I)*W(L)*Q(I,J)*Q(I,J)   
      T2= 3.*W(I)*W(L)*P(I)*P(J)+3.*W2(J)*P(I)*Q(I,L)+6.*W(J)*W(L)*P(I  
     & )*Q(I,J)-7.*P(I)*Q(I,J)*Q(J,L)-2.*P(I)*Q(I,L)*P(J)-6.*Q(I,J)*Q(  
     & I,J)*Q(I,L)                                                      
      H5=T1+T2                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H6(I,J,L)                                                
C                                                                       
C     HERMIT POLYNOM OF SIXTH ORDER      H6 = H(I,I,J,J,L,L)            
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      T1= W2(I)*W2(J)*W2(L)-W2(I)*W2(J)*P(L)-W2(I)*W2(L)*P(J)-W2(J)*W2  
     & (L)*P(I)-4.*W2(I)*W(J)*W(L)*Q(J,L)-4.*W(I)*W2(J)*W(L)*Q(I,L)-4.  
     & *W(I)*W(J)*W2(L)*Q(I,J)+W2(I)*P(J)*P(L)+W2(J)*P(I)*P(L)+W2(L)*P  
     & (I)*P(J)+2.*W2(I)*Q(J,L)*Q(J,L)+2.*W2(J)*Q(I,L)*Q(I,L)           
      T2=2.*W2(L)*Q(I,J)*Q(I,J)+4.*W(I)*W(J)*Q(I,J)*P(L)+4.*W(I)*W(L)*  
     & Q(I,L)*P(J)+4.*W(J)*W(L)*Q(J,L)*P(I)+8.*W(I)*W(J)*Q(I,L)*Q(J,L)  
     & +8.*W(I)*W(L)*Q(I,J)*Q(J,L)+8.*W(J)*W(L)*Q(I,L)*Q(I,J)           
      T1=T1+T2                                                          
      T2=-P(I)*P(J)*P(L)-2.*P(I)*Q(J,L)*Q(J,L)-2.*P(J)*                 
     & Q(I,L)*Q(I,L)-2.*P(L)*Q(I,J)*Q(I,J)-8.*Q(I,J)*Q(I,L)*Q(J,L)      
      H6=T1+T2                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION H7(I,J)                                                  
C                                                                       
C     HERMIT POLYNOM OF SIXTH RANK      H7 = H(I,I,I,J,J,J)             
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      H7= W3(I)*W3(J)-3.*W3(I)*W(J)*P(J)-9.*W2(I)*W2(J)*Q(I,J)-3.*W(I)  
     & *W3(J)*P(I)+9.*W2(I)*Q(I,J)*P(J)+18.*W(I)*W(J)*Q(I,J)*Q(I,J)+9.  
     & *W(I)*W(J)*P(I)*P(J)+9.*W2(J)*P(I)*Q(I,J)-9.*Q(I,J)*P(I)*P(J)-6  
     & .*Q(I,J)**3                                                      
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE SHERE                                                  
C                                                                       
C---- CALCULATION OF OUTPUT ANGLE OF THE POINT ARRAY.                   
C                                                                       
      COMMON /PRINT/WINKEL(3),ISELEC,ISECT,IS,IHILF,LSPCH,GAM           
C                                                                       
      GAM=WINKEL(3)                                                     
      IF(ISELEC.GE.2)GAM=FLOAT(ISELEC)                                  
      IF(ISELEC.GE.2)GO TO 113                                          
      IF(ISECT.EQ.1.OR.ISECT.EQ.5)GAM=WINKEL(2)                         
      IF(ISECT.EQ.3.OR.ISECT.EQ.4)GAM=WINKEL(1)                         
  113 IF(GAM.LT.30..OR.GAM.GT.160)GAM=90.                               
      IF (GAM-90.0) 114,115,116                                         
  115 IS=0                                                              
      GO TO 117                                                         
  116 HILF=GAM-90.0                                                     
      IS=1                                                              
  118 IF(LSPCH.NE.0)DEL=30.0                                            
      IF(LSPCH.EQ.0)DEL=15.0                                            
      HILF=HILF/DEL                                                     
      IHILF=INT(HILF)                                                   
      DIF=HILF-FLOAT(IHILF)                                             
      IF (DIF.GE.0.5) IHILF=IHILF+1                                     
      IF(IHILF.EQ.0)IS=0                                                
      GOTO 117                                                          
  114 HILF=ABS(GAM-90)                                                  
      IS=-1                                                             
      GO TO 118                                                         
  117 RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ARB(ISECT,CELL,ALF,BET,GAM,AG,BG,CG,JACHS)             
C                                                                       
C---- CHANGE SEQUENCE OF LATTICE PARAMETERS CORRESPONDING TO            
C     SELECTED SECTIONS.                                                
C                                                                       
      DIMENSION CELL(6)                                                 
      IF(ISECT)1,1,30                                                   
   30 GOTO(2,3,4,5,6),ISECT                                             
    1 CELL(1)=AG                                                        
      CELL(2)=BG                                                        
      CELL(3)=CG                                                        
      CELL(4)=ALF                                                       
      CELL(5)=BET                                                       
      CELL(6)=GAM                                                       
      JACHS=0                                                           
      RETURN                                                            
    2 CELL(1)=AG                                                        
      CELL(2)=CG                                                        
      CELL(3)=BG                                                        
      CELL(4)=ALF                                                       
      CELL(5)=GAM                                                       
      CELL(6)=BET                                                       
      JACHS=1                                                           
      RETURN                                                            
    3 CELL(1)=BG                                                        
      CELL(2)=AG                                                        
      CELL(3)=CG                                                        
      CELL(4)=BET                                                       
      CELL(5)=ALF                                                       
      CELL(6)=GAM                                                       
      JACHS=2                                                           
      RETURN                                                            
    4 CELL(1)=BG                                                        
      CELL(2)=CG                                                        
      CELL(3)=AG                                                        
      CELL(4)=BET                                                       
      CELL(5)=GAM                                                       
      CELL(6)=ALF                                                       
      JACHS=3                                                           
      RETURN                                                            
    5 CELL(1)=CG                                                        
      CELL(2)=BG                                                        
      CELL(3)=AG                                                        
      CELL(4)=GAM                                                       
      CELL(5)=BET                                                       
      CELL(6)=ALF                                                       
      JACHS=4                                                           
      RETURN                                                            
    6 CELL(1)=CG                                                        
      CELL(2)=AG                                                        
      CELL(3)=BG                                                        
      CELL(4)=GAM                                                       
      CELL(5)=ALF                                                       
      CELL(6)=BET                                                       
      JACHS=5                                                           
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION HJK(I)                                                   
C                                                                       
C     CALCULATION OF SECOND ORDER HERMIT POLYNOMS IN THREE DIM.         
C                                                                       
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)         
      COMMON /INVBET/ Q(3,3)                                            
C                                                                       
      IF(I.GT.3) GOTO 1                                                 
      HJK =W(I)*W(I)-P(I)                                               
      RETURN                                                            
    1 IF(I.GT.5) GOTO 2                                                 
      J=I-2                                                             
      HJK=W(1)*W(J)-Q(1,J)                                              
      RETURN                                                            
    2 HJK=W(2)*W(3)-Q(2,3)                                              
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE TRANSF(ISECT,I1,I2,I3,M3)                              
C                                                                       
C                                                                       
C     CELL TRANSFORMATION IN ACCORDANCE TO SELECTED SECTIONS            
C                                                                       
      REAL M3                                                           
C                                                                       
      GOTO (2,3,4,5,6),ISECT                                            
C                                                                       
C     XZY DOWN                                                          
C                                                                       
   2  I2=3                                                              
      I3=2                                                              
      M3=1.                                                             
      RETURN                                                            
C                                                                       
C     YXZ                                                               
C                                                                       
   3  I1=2                                                              
      I2=1                                                              
      M3=1.                                                             
      RETURN                                                            
C                                                                       
C     YZX                                                               
C                                                                       
   4  I1=2                                                              
      I2=3                                                              
      I3=1                                                              
      RETURN                                                            
C                                                                       
C     ZYX                                                               
C                                                                       
   5  I1=3                                                              
      I3=1                                                              
      M3=1.                                                             
      RETURN                                                            
C                                                                       
C     ZXY                                                               
C                                                                       
   6  I1=3                                                              
      I2=1                                                              
      I3=2                                                              
      RETURN                                                            
C                                                                       
      END                                                               
C                                                                       
      FUNCTION RANK5(ER,JFIVE,II)                                       
C                                                                       
C     CALCULATION OF                                                    
C                                                                       
C                I J K L M                                              
C     RANK5 = ER           * HERMIT5                                    
C                                    I J K L M                          
C                                                                       
C                                                                       
C REMARK=  THE  FOLLOWING PROGRAMMING SEEMS TO BE  VERY  COMPLICATED AND
C          INEFFECTIVE.   BUT I SUPOSSE NOBODY WILL CALCULATED FIFTH  OR
C          SIXTH CUMULANTS  FOR ATOMS IN GENERAL POSITIONS. THAT   MEANS
C          THAT IN MOST OF ALL CASES A LOT OF CUMULANTS HAVE TO BE ZERO.
C          BECAUSE OF COMPLEXITY OF HERMIT POLYNOMS OF  FIFTH OR  HIGHER
C          ORDER, THIS WAY  OF PROGRAMMING SAVES A REASONABLE AMOUNT  OF
C          PROCESSING TIME.     TIME IS MONEY||                         
C                                                                       
       DIMENSION ER(21,6),JFIVE(21)                                     
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)
      COMMON /INVBET/ Q(3,3)
C                                                                       
C      INITIALIZE                                                       
C                                                                       
       CALL START                                                       
       CUM5 = 0.0                                                       
C                                                                       
       IF(JFIVE(1))2900,2901,2900                                       
 2900  CUM5 =        ER(1,II)*HH1(1)                                    
 2901  IF(JFIVE(2))2902,2903,2902                                       
 2902  CUM5 = CUM5 + ER(2,II)*HH1(2)                                    
 2903  IF(JFIVE(3))2904,2905,2904                                       
 2904  CUM5 = CUM5 + ER(3,II)*HH1(3)                                    
 2905  TEMP = 0.0                                                       
       IF(JFIVE(4))2906,2907,2906                                       
 2906  TEMP =        ER(4,II)*HH2(1,2)                                  
 2907  IF(JFIVE(5))2908,2909,2908                                       
C                                                                       
 2908  TEMP = TEMP + ER(5,II)*HH2(1,3)                                  
 2909  IF(JFIVE(13))2910,2911,2910                                      
 2910  TEMP = TEMP + ER(13,II)*HH2(2,1)                                 
 2911  IF(JFIVE(18))2912,2913,2912                                      
 2912  TEMP = TEMP + ER(18,II)*HH2(2,3)                                 
 2913  IF(JFIVE(17))2914,2915,2914                                      
 2914  TEMP = TEMP + ER(17,II)*HH2(3,1)                                 
 2915  IF(JFIVE(21))2916,2917,2916                                      
 2916  TEMP = TEMP + ER(21,II)*HH2(3,2)                                 
 2917  CUM5 = CUM5 + TEMP*5.                                            
       TEMP = 0.0                                                       
       IF(JFIVE(6))2918,2919,2918                                       
 2918  TEMP =        ER(6,II) * HH3(1,2)                                
 2919  IF(JFIVE(8))2920,2921,2920                                       
 2920  TEMP = TEMP + ER(8,II) * HH3(1,3)                                
 2921  IF(JFIVE(19))2922,2923,2922                                      
 2922  TEMP = TEMP + ER(19,II) * HH3(2,3)                               
 2923  IF(JFIVE(9)) 2924,2925,2924                                      
 2924  TEMP = TEMP + ER( 9,II) * HH3(2,1)                               
 2925  IF(JFIVE(12))2926,2927,2926                                      
 2926  TEMP = TEMP + ER(12,II) * HH3(3,1)                               
 2927  IF(JFIVE(20))2928,2929,2928                                      
 2928  TEMP = TEMP + ER(20,II) * HH3(3,2)                               
 2929  CUM5 = CUM5 +TEMP*10.                                            
       TEMP = 0.0                                                       
       IF(JFIVE( 7))2930,2931,2930                                      
 2930  TEMP =        ER( 7,II) * HH4(1,2,3)                             
 2931  IF(JFIVE(14))2932,2933,2932                                      
 2932  TEMP = TEMP + ER(14,II) * HH4(2,1,3)                             
 2933  IF(JFIVE(16))2934,2935,2934                                      
 2934  TEMP = TEMP + ER(16,II) * HH4(3,1,2)                             
 2935  CUM5 = CUM5 + TEMP*20.                                           
       TEMP = 0.0                                                       
       IF(JFIVE(10))2936,2937,2936                                      
 2936  TEMP =        ER(10,II) * HH5(1,2,3)                             
 2937  IF(JFIVE(11))2938,2939,2938                                      
 2938  TEMP = TEMP + ER(11,II) * HH5(1,3,2)                             
 2939  IF(JFIVE(15))2940,2941,2940                                      
 2940  TEMP = TEMP + ER(15,II) * HH5(2,3,1)                             
 2941  CUM5 = CUM5 + TEMP*30.                                           
C                                                                       
       RANK5= CUM5 * 1.E-5                                              
C                                                                       
       RETURN                                                           
       END                                                              
C                                                                       
       FUNCTION RANK6(FR,JSIX,NO5,II)                                   
C                                                                       
C      CALCULATE                                                        
C                                                                       
C                I J K L M N                                            
C      RANK6 = FR            * HERMIT6                                  
C                                      I J K L M N                      
C                                                                       
C                                                                       
C                                                                       
       DIMENSION FR(28,6),JSIX(28)                                      
      COMMON /WQ/ W(3),W2(3),W3(3),W4(3),W5(3),P(3),P2(3),P3(3)
      COMMON /INVBET/ Q(3,3)
C                                                                       
C----  INITIALIZE, IF NOT DONE BEFORE.                                  
C                                                                       
       IF(NO5.EQ.0) CALL START                                          
             CUM6 = 0.                                                  
C                                                                       
C----  IN CASE OF ASTONISHMENT SEE REMARK IN FUNCTION RANK5.            
C                                                                       
       DO 2011 IDD=1,3                                                   
      ID = IDD
       IF(JSIX(ID).EQ.0) GOTO 2011                                      
       CUM6 = CUM6 + FR(ID,II)*H1(ID)                                   
 2011  CONTINUE                                                         
       TEMP = 0.0                                                       
       IF(JSIX( 4))2012,2013,2012                                       
 2012  TEMP =        FR(4,II)*H2(1,2)                                   
 2013  IF(JSIX( 5))2014,2015,2014                                       
 2014  TEMP = TEMP + FR(5,II)*H2(1,3)                                   
 2015  IF(JSIX(18))3016,3017,3016                                       
 3016  TEMP = TEMP + FR(18,II)*H2(2,1)                                  
 3017  IF(JSIX(24))2016,2017,2016                                       
 2016  TEMP = TEMP + FR(24,II)*H2(2,3)                                  
 2017  IF(JSIX(23))2018,2019,2018                                       
 2018  TEMP = TEMP + FR(23,II)*H2(3,1)                                  
 2019  IF(JSIX(28))2020,2021,2020                                       
 2020  TEMP = TEMP + FR(28,II)*H2(3,2)                                  
 2021  CUM6 = CUM6 + TEMP*6.                                            
       TEMP = 0.0                                                       
       IF(JSIX( 6))2022,2023,2022                                       
 2022  TEMP =        FR( 6,II)*H3(1,2)                                  
 2023  IF(JSIX( 8))2024,2025,2024                                       
 2024  TEMP = TEMP + FR( 8,II)*H3(1,3)                                  
 2025  IF(JSIX(25))2026,2027,2026                                       
 2026  TEMP = TEMP + FR(25,II)*H3(2,3)                                  
 2027  IF(JSIX(13))2028,2029,2028                                       
 2028  TEMP = TEMP + FR(13,II)*H3(2,1)                                  
 2029  IF(JSIX(17))2030,2031,2030                                       
 2030  TEMP = TEMP + FR(17,II)*H3(3,1)                                  
 2031  IF(JSIX(27))2032,2033,2032                                       
 2032  TEMP = TEMP + FR(27,II)*H3(3,2)                                  
 2033  CUM6 = CUM6 + TEMP*15.                                           
       TEMP = 0.0                                                       
       IF(JSIX( 7))2034,2035,2034                                       
 2034  TEMP =        FR( 7,II)*H4(1,2,3)                                
 2035  IF(JSIX(19))2036,2037,2036                                       
 2036  TEMP = TEMP + FR(19,II)*H4(2,1,3)                                
 2037  IF(JSIX(22)) 2038,2039,2038                                      
 2038  TEMP = TEMP + FR(22,II)*H4(3,1,2)                                
 2039  CUM6 =CUM6 + TEMP*30.                                            
       TEMP =0.0                                                        
       IF(JSIX( 9))2040,2041,2040                                       
 2040  TEMP =        FR( 9,II)*H7(1,2)                                  
 2041  IF(JSIX(12))2042,2043,2042                                       
 2042  TEMP = TEMP + FR(12,II)*H7(1,3)                                  
 2043  IF(JSIX(26))2044,2045,2044                                       
 2044  TEMP = TEMP + FR(26,II)*H7(2,3)                                  
 2045  CUM6 = CUM6 + TEMP*20.                                           
       TEMP = 0.0                                                       
       IF(JSIX(10))2046,2047,2046                                       
 2046  TEMP =        FR(10,II)*H5(1,2,3)                                
 2047  IF(JSIX(11))2048,2049,2048                                       
 2048  TEMP = TEMP + FR(11,II)*H5(1,3,2)                                
 2049  IF(JSIX(14))2050,2051,2050                                       
 2050  TEMP = TEMP + FR(14,II)*H5(2,1,3)                                
 2051  IF(JSIX(16))2052,2053,2052                                       
 2052  TEMP = TEMP + FR(16,II)*H5(3,1,2)                                
 2053  IF(JSIX(20))2054,2055,2054                                       
 2054  TEMP = TEMP + FR(20,II)*H5(2,3,1)                                
 2055  IF(JSIX(21))2056,2057,2056                                       
 2056  TEMP = TEMP + FR(21,II)*H5(3,2,1)                                
 2057  CUM6 = CUM6 + TEMP*60.                                           
       IF(JSIX(15))2058,2059,2058                                       
 2058  CUM6 = CUM6 + FR(15,II)*H6(1,2,3)*90.                            
 2059  RANK6= CUM6*1.E-6                                                
C                                                                       
       RETURN                                                           
       END                                                              
                          SUBROUTINE JOHN                               
     &  (W,Q,ITHREE,JTHREE,CR,IFOUR,JFOUR,DR,ALF,II,CUM3,CUM4,NORM)     
C                                                                       
C     THIS ROUTINE CALCULATES THE FOURIER TRANSFORM OF THE  STRUCTURE   
C     FACTOR  BASED ON THE  ALPHA FORMALISM.  THIS  STRUCTURE  FACTOR   
C     EQUATION IS IN THE PROGRAMME SYSTEM PROMETHEUS CALLED  MODIFIED   
C     GRAM-CHARLIER SERIES EXPANSION. (MODEL NO. FOUR)                  
C     THE FOURIER TRANSFORM HAD BEEN PUBLISHED BY JOHNSON IN  THE OAK   
C     RIDGE REPORT 1980  DOE/TIC-11068 THEME=THERMAL MOTION ANALYSIS.   
C     THE EQUATION (9) OF THIS REPORT IS INCORRECTLY NORMALIZED. THIS   
C     WILL BE DONE HERE IN ADDITION. THE FOURIER TRANSFORM ARE RECAL-   
C     CULATED BY THE AUTHOR OF THIS PROGRAMME AND FIND TO BE CORRECT.   
C                                                                       
C                    UDO ZUCKER  NOVEMBER 1980                          
C                                                                       
C     IN THE ROUTINE 3 QUANTITYS WILL BE CALCULATED. CUM3,CUM4,NORM.    
C                                                                       
      REAL NORM                                                         
C                                                                       
      DIMENSION W(3),Q(3,3),CR(10,60),DR(15,60),JTHREE(10),JFOUR(15)    
C                                                                       
      W1    =  W(1)                                                     
      W2    =  W(2)                                                     
      W3    =  W(3)                                                     
      W11   = W1*W1                                                     
      W22   = W2*W2                                                     
      W33   = W3*W3                                                     
C                                                                       
C---- NORM IS A REAL VAIABLE WHICH CARRIES INFO. ABOUT NORMALIZATION    
C                                                                       
C---- THE VARIABLE ALPHA OF THE PROGRAMME SYSTEM PROMETHEUS IS ONE      
C     HALF OF THE VARIABLE ALPHA BY JOHNSON. HERE A WILL BE SET TO      
C     TWO ALF, TO USE THE FORMULAS BY JOHNSON. (Y  HAS BEEN REPLA-      
C     CED BY  W  TO BE CONSISTENT WITH INTERNATIONAL TABLES OF X -      
C     RAY CRYSTALLOGRAPHY NO.4)                                         
C                                                                       
      A     = 2.*ALF                                                    
      OMA   = 1.-A                                                      
      OMAQ2 = OMA**2                                                    
C                                                                       
C---- NO THIRD ORDER PRODUCT,IF ITHREE EQUAL ZERO.                      
C                                                                       
      IF(ITHREE.EQ.0) GOTO 3000                                         
C                                                                       
C---- INITIALIZE                                                        
C                                                                       
      TERM  = 0.                                                        
      TERMZ = 0.                                                        
C                                                                       
      IF(JTHREE(1))1,2,1                                                
    1 TERM1 = W11*W1                                                    
      TERM2 = 3.*W1*Q(1,1)*OMA                                          
      TERM  = (TERM1-TERM2)*CR(1,II)                                    
    2 IF(JTHREE(2))3,4,3                                                
    3 TERM1 = W22*W2                                                    
      TERM2 = 3.*W2*Q(2,2)*OMA                                          
      TERM  = TERM + (TERM1-TERM2)*CR(2,II)                             
    4 IF(JTHREE(3))5,6,5                                                
    5 TERM1 = W33*W3                                                    
      TERM2 = 3.*W3*Q(3,3)*OMA                                          
      TERM  = TERM + (TERM1-TERM2)*CR(3,II)                             
    6 IF(JTHREE(4))7,8,7                                                
    7 TERM1 = W11*W2                                                    
      TERM2 = (2.*W1*Q(1,2)+W2*Q(1,1))*OMA                              
      TERMZ = (TERM1-TERM2)*CR(4,II)                                    
    8 IF(JTHREE(5))9,10,9                                               
    9 TERM1 = W22*W1                                                    
      TERM2 = (W1*Q(2,2)+2.*W2*Q(1,2))*OMA                              
      TERMZ = TERMZ + (TERM1-TERM2)*CR(5,II)                            
   10 IF(JTHREE(6))11,12,11                                             
   11 TERM1 = W11*W3                                                    
      TERM2 = (2.*W1*Q(1,3)+W3*Q(1,1))*OMA                              
      TERMZ = TERMZ + (TERM1-TERM2)*CR(6,II)                            
   12 IF(JTHREE(7))13,14,13                                             
   13 TERM1 = W1*W33                                                    
      TERM2 = (W1*Q(3,3)+2.*W3*Q(1,3))*OMA                              
      TERMZ = TERMZ + (TERM1-TERM2)*CR(7,II)                            
   14 IF(JTHREE(8))15,16,15                                             
   15 TERM1 = W22*W3                                                    
      TERM2 = (2.*W2*Q(2,3)+W3*Q(2,2))*OMA                              
      TERMZ = TERMZ + (TERM1-TERM2)*CR(8,II)                            
   16 IF(JTHREE(9))17,18,17                                             
   17 TERM1 = W2*W33                                                    
      TERM2 = (W2*Q(3,3)+2.*W3*Q(2,3))*OMA                              
      TERMZ = TERMZ + (TERM1-TERM2)*CR(9,II)                            
   18 TERM  = TERM + TERMZ*3.0                                          
      TERMZ = 0.0                                                       
      IF(JTHREE(10))19,20,19                                            
   19 TERM1 = W1*W2*W3                                                  
      TERM2 = (W1*Q(2,3)+W2*Q(1,3)+W3*Q(1,2))*OMA                       
      TERMZ = (TERM1-TERM2)*CR(10,II)*6.0                               
   20 CUM3  = TERM + TERMZ                                              
      CUM3  = CUM3 * 1.E-3                                              
C                                                                       
C---- IF IFOUR=0 FOURTH RANK PRODUCTS ARE NOT REQUIRED.                 
C                                                                       
 3000 IF(IFOUR.EQ.0) GOTO 4000                                          
C                                                                       
      TERM  = 0.                                                        
      TERMZ = 0.                                                        
C                                                                       
      IF(JFOUR(1))101,102,101                                           
  101 TERM1 = W11*W11                                                   
      TERM2 = 6.*W11*Q(1,1)*OMA                                         
      TERM3 = 3.*Q(1,1)**2                                              
      NORM  = TERM3                                                     
      TERM3 = TERM3*OMAQ2                                               
      TERM  = (TERM1-TERM2+TERM3)*DR(1,II)                              
  102 IF(JFOUR(2))103,104,103                                           
  103 TERM1 = W22*W22                                                   
      TERM2 = 6.*W22*Q(2,2)*OMA                                         
      TERM3 = 3.*Q(2,2)**2                                              
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERM  = TERM + (TERM1-TERM2+TERM3)*DR(2,II)                       
  104 IF(JFOUR(3))105,106,105                                           
  105 TERM1 = W33*W33                                                   
      TERM2 = 6.*W33*Q(3,3)*OMA                                         
      TERM3 = 3.*Q(3,3)**2                                              
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERM  = TERM + (TERM1-TERM2+TERM3)*DR(3,II)                       
  106 IF(JFOUR(4))107,108,107                                           
  107 TERM1 = W11*W1*W2                                                 
      TERM2 = (W11*Q(1,2)*3.+W1*W2*Q(1,1)*3.)*OMA                       
      TERM3 = 3.*Q(1,1)*Q(1,2)                                          
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = (TERM1-TERM2+TERM3)*DR(4,II)                              
  108 IF(JFOUR(5))109,110,109                                           
  109 TERM1 = W11*W1*W3                                                 
      TERM2 = (W11*Q(1,3)*3.+W1*W3*Q(1,1)*3.)*OMA                       
      TERM3 = 3.*Q(1,1)*Q(1,3)                                          
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(5,II)                      
  110 IF(JFOUR(9))111,112,111                                           
  111 TERM1 = W22*W2*W1                                                 
      TERM2 = (W22*Q(1,2)*3.+W2*W1*Q(2,2)*3.)*OMA                       
      TERM3 = 3.*Q(2,2)*Q(1,2)                                          
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(9,II)                      
  112 IF(JFOUR(12))113,114,113                                          
  113 TERM1 = W33*W3*W1                                                 
      TERM2 = (W33*Q(1,3)*3.+W1*W3*Q(3,3)*3.)*OMA                       
      TERM3 = 3.*Q(3,3)*Q(1,3)                                          
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(12,II)                     
  114 IF(JFOUR(13))115,116,115                                          
  115 TERM1 = W22*W2*W3                                                 
      TERM2 = (W22*Q(2,3)*3.+W2*W3*Q(2,2)*3.)*OMA                       
      TERM3 = 3*Q(2,2)*Q(2,3)                                           
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(13,II)                     
  116 IF(JFOUR(15))117,118,117                                          
  117 TERM1 = W33*W3*W2                                                 
      TERM2 = (W33*Q(2,3)*3.+W2*W3*Q(3,3)*3.)*OMA                       
      TERM3 = 3.*Q(3,3)*Q(2,3)                                          
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(15,II)                     
  118 TERM  = TERM + TERMZ*4.0                                          
      TERMZ = 0.0                                                       
      IF(JFOUR(6))119,120,119                                           
  119 TERM1 = W11*W22                                                   
      TERM2 = (W11*Q(2,2)+4.*W1*W2*Q(1,2)+W22*Q(1,1))*OMA               
      TERM3 = Q(1,1)*Q(2,2)+2.*Q(1,2)**2                                
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = (TERM1-TERM2+TERM3)*DR(6,II)                              
  120 IF(JFOUR(8))121,122,121                                           
  121 TERM1 = W11*W33                                                   
      TERM2 = (W11*Q(3,3)+4.*W1*W3*Q(1,3)+W33*Q(1,1))*OMA               
      TERM3 = Q(1,1)*Q(3,3)+2.*Q(1,3)**2                                
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(8,II)                      
  122 IF(JFOUR(14))123,124,123                                          
  123 TERM1 = W33*W22                                                   
      TERM2 = (W33*Q(2,2)+4.*W3*W2*Q(2,3)+W22*Q(3,3))*OMA               
      TERM3 = Q(3,3)*Q(2,2)+2.*Q(2,3)**2                                
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(14,II)                     
  124 TERM  = TERM + TERMZ*6.0                                          
      TERMZ =0.0                                                        
      IF(JFOUR(10))125,126,125                                          
  125 TERM1 = W22*W3*W1                                                 
      TERM2 = (W22*Q(1,3)+2.*W2*W3*Q(1,2)+2.*W2*W1*Q(2,3) +             
     &        W3*W1*Q(2,2)) *OMA                                        
      TERM3 = Q(2,2)*Q(1,3)+2.*Q(2,3)*Q(1,2)                            
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = (TERM1-TERM2+TERM3)*DR(10,II)                             
  126 IF(JFOUR(11))127,128,127                                          
  127 TERM1 = W1*W2*W33                                                 
      TERM2 = (W1*W2*Q(3,3)+W2*W3*Q(1,3)*2.+W33*Q(1,2) +                
     &        2.*W1*W3*Q(2,3)) *OMA                                     
      TERM3 = Q(1,2)*Q(3,3)+2.*Q(1,3)*Q(2,3)                            
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(11,II)                     
  128 IF(JFOUR(7))129,130,129                                           
  129 TERM1 = W11*W2*W3                                                 
      TERM2 = (W11*Q(2,3)+2.*W1*W2*Q(1,3)+2.*W1*W3*Q(1,2) +             
     &        W2*W3*Q(1,1)) *OMA                                        
      TERM3 = Q(1,1)*Q(2,3)+2.*Q(1,2)*Q(1,3)                            
      NORM  = NORM + TERM3                                              
      TERM3 = TERM3*OMAQ2                                               
      TERMZ = TERMZ + (TERM1-TERM2+TERM3)*DR(7,II)                      
  130 CUM4  = TERM + TERMZ*12.0                                         
      CUM4  = CUM4 * 1.E-4                                              
C                                                                       
C                                                                       
      NORM = NORM * A*A                                                 
C                                                                       
C                                                                       
 4000 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C                                                                       
                         SUBROUTINE PLANE                               
     &(P1,P2,P3,IX,IY,IXI,IXF,IYI,IYF,ITYPE,G,X,Y,Z,AGA,BGB,A,B,IERR)   
C                                                                       
C     THIS  PROGRAMME CALCULATES THE UNIT VECTORS OF A SPECIAL  PLANE   
C           DEFINED BY THREE POINTS AND A GRID COORDINATES.             
C                                                                       
C                   BY   ROLF BACHMANN   DEC. 1980                      
C                                                                       
C                                                                       
C                         STUTTGART  FRG                                
C                                                                       
C                                                                       
C----INPUT    THREE  POINTS P1,P2,P3, THE GRID IX*IY                    
C             LIMITS OF THE PLANE IXI, IXF, IYI, IYF                    
C             TYPE ITYPE, AND THE METRIC G                              
C----OUTPUT   ORIGIN  FOR THE PLANE, X, Y, Z, LENGTHS AGA, BGB          
C             UNIT VECTORS SPANNING THE PLANE A, B                      
C             IERR, IF THE THREE POINTS LIE ON A STRAIGHT LINE          
C                                                                       
      DIMENSION P1(3),P2(3),P3(3),A(3),B(3),G(3,3)                      
C                                                                       
      COMMON /IWRITE/AL,BL,CL,CAL,CBE,CGA,AM(3),C,CD,ALANG,BLANG,IL,ILY 
      COMMON /APLANE/AA(3)                                              
C                                                                       
C---- CALCULATING CELL CONSTANTS FROM G                                 
C                                                                       
      AL=SQRT(G(1,1))                                                   
      BL=SQRT(G(2,2))                                                   
      CL=SQRT(G(3,3))                                                   
      CAL=ACOS(G(2,3)/BL/CL)*57.295779                                  
      CBE=ACOS(G(1,3)/AL/CL)*57.295779                                  
      CGA=ACOS(G(1,2)/AL/BL)*57.295779                                  
C                                                                       
C---- VECTOR A=P2-P1, VECTOR B=P3-P1                                    
C                                                                       
      A(1)=P2(1)-P1(1)                                                  
      A(2)=P2(2)-P1(2)                                                  
      A(3)=P2(3)-P1(3)                                                  
      B(1)=P3(1)-P1(1)                                                  
      B(2)=P3(2)-P1(2)                                                  
      B(3)=P3(3)-P1(3)                                                  
C                                                                       
C---- PRODUCT OF VECTORS A*A, A*B, B*B                                  
C                                                                       
      AGB=A(1)*(G(1,1)*B(1)+G(1,2)*B(2)+G(1,3)*B(3))+                   
     &    A(2)*(G(2,1)*B(1)+G(2,2)*B(2)+G(2,3)*B(3))+                   
     &    A(3)*(G(3,1)*B(1)+G(3,2)*B(2)+G(3,3)*B(3))                    
      AGA=A(1)*(G(1,1)*A(1)+2*G(1,2)*A(2)+2*G(1,3)*A(3))+               
     &    A(2)*(G(2,2)*A(2)+2*G(2,3)*A(3))+A(3)*G(3,3)*A(3)             
      BGB=B(1)*(G(1,1)*B(1)+2*G(1,2)*B(2)+2*G(1,3)*B(3))+               
     &    B(2)*(G(2,2)*B(2)+2*G(2,3)*B(3))+B(3)*G(3,3)*B(3)             
      ALANG=SQRT(AGA)                                                   
C                                                                       
C---- ARE THE POINTS LINEAR DEPENDENT?                                  
C                                                                       
      IF((AGA*BGB).LT.1.E-5)GOTO 149                                    
C                                                                       
C---- CALCULATING THE ANGLE BETWEEN A AND B, IF THE COSINE              
C---- IS VERY CLOSE TO OR 1 -1 IERR=1                                   
C                                                                       
      CD= AGB/SQRT(AGA*BGB)                                             
      C=ABS(CD)-1.                                                      
      IF(C.LT.-1.E-5)GOTO 150                                           
  149 IERR=1                                                            
      RETURN                                                            
C                                                                       
C---- THE VECTOR PERPENDICULAR TO P1-P2 THROUGH P3 IS CALCULATED        
C                                                                       
  150 C=AGB/AGA                                                         
      AM(1)=B(1)-A(1)*C                                                 
      AM(2)=B(2)-A(2)*C                                                 
      AM(3)=B(3)-A(3)*C                                                 
C                                                                       
C---- IF P3 IS CENTERPOINT DOUBLE AM, AM IS THE Y AXIS OF THE PLANE     
C                                                                       
      IF(ITYPE.NE.1)GOTO 200                                            
      AM(1)=AM(1)*2                                                     
      AM(2)=AM(2)*2                                                     
      AM(3)=AM(3)*2                                                     
C                                                                       
C---- CALCULATING UNIT VECTORS                                          
C                                                                       
  200 AA(1)=AM(1)                                                       
      AA(2)=AM(2)                                                       
      AA(3)=AM(3)                                                       
      AIX=FLOAT(IX)                                                     
      AIY=FLOAT(IY)                                                     
      IL=IX                                                             
      ILY=IY                                                            
      A(1)=A(1)/AIX                                                     
      A(2)=A(2)/AIX                                                     
      A(3)=A(3)/AIX                                                     
      B(1)=AM(1)                                                        
      B(2)=AM(2)                                                        
      B(3)=AM(3)                                                        
      BGB=B(1)*(G(1,1)*B(1)+2*G(1,2)*B(2)+2*G(1,3)*B(3))+               
     &    B(2)*(G(2,2)*B(2)+2*G(2,3)*B(3))+B(3)*G(3,3)*B(3)             
      BLANG=SQRT(BGB)                                                   
      B(1)=B(1)/AIY                                                     
      B(2)=B(2)/AIY                                                     
      B(3)=B(3)/AIY                                                     
C                                                                       
C---- ORIGIN OF THE PLANE                                               
C                                                                       
      AIXI=FLOAT(IXI)                                                   
      AIYI=FLOAT(IYI)                                                   
      X=P1(1)+A(1)*AIXI+B(1)*AIYI                                       
      Y=P1(2)+A(2)*AIXI+B(2)*AIYI                                       
      Z=P1(3)+A(3)*AIXI+B(3)*AIYI                                       
      AGA=SQRT(AGA)                                                     
      BGB=SQRT(BGB)                                                     
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
                          SUBROUTINE SRBPLA                             
     &    (G,P1,P2,P3,X,Y,Z,A,B,IOUT,IERR,ITYPE,MP,PP,IX,IY,IA)         
C                                                                       
C                                                                       
C     THIS  PROGRAMM  WRITES THE  CONSTANTS  OF THE PLANE, CALCULATED   
C     BY THE SUROUTINE PLANE, AS WELL AS THE CELL CONSTANTS AND G. IT   
C     DOES ALSO GIVE YOU A ROUGH SKETCH OF HOW P1, P2 AND P3  LIE  IN   
C     THIS PLANE AND  WRITES THEIR  FRACTIONAL  COORDINATES,  IT ALSO   
C     SHOWS YOU THE POSITION OF AN ATOM.                                
C                                                                       
C                  BY ROLF BACHMANN  DECEMBER 1980                      
C                                                                       
C     INPUT      THE  ORIGIN OF THE PLANE  X,Y,Z,THE POINTS P1,P2,P3    
C                THE  METRIC G,  VECTORS SPANNING THE PLANE A,B, THE    
C                OUTPUTUNIT IOUT, IERR IF P1, P2,P3 LIE ON A STRAIGT    
C                LINE. ITYPE SPECIFIES THE TYPE OF THE DESIRED PLANE    
C                THE GRID IX*IY,  PP THE GRID COORDINATES OF AN ATOM    
C                AND  IA IF YOU WANT  THE ATOM TO BE DRAWN.  ALNU IS    
C                JUST AN OUTPUT ARRAY.                                  
C                                                                       
C                                                                       
C---- ALNU IS EQUIVALENT TO ANOTHER ARRAY NOT USED AT THE MOMENT.       
C                                                                       
      DOUBLE PRECISION NAME                                             
C                                                                       
      DIMENSION ALNU(80,22),P1(3),P2(3),P3(3),G(3,3),A(3),B(3)          
      DIMENSION PP1(2),PP2(2),PP3(2),PP(3),NAME(2),ARRAY(10000)
      Dimension qs(10000),iplus(10000),iminus(10000)                   
      dimension qp(10000),qm(10000)   
C                                                                       
      COMMON /IWRITE/AL,BL,CL,CAL,CBE,CGA,AM(3),C,CD,ALANG,BLANG,IL,ILY 
      COMMON /BACH/ARRAY,qp,qm,qs,iplus,iminus                                             
C                                                                       
C    EQUIVALENCE (ARRAY,ALNU)                                          
C                                                                       
      DATA AI,AB,AS,AP,AMI,AP1,AP2 /1HI,1H ,1H*,1HP,1H-,1H1,1H2/        
      DATA ASS,AYY,AXX,AA,AP3 /1HS,1HY,1HX,1HA,1H3/                     
      DATA AMM,ATT,AOO,NAME /1HM,1HT,1HO,8H CORNER ,8H CENTER /         
C                                                                       
   11 FORMAT(1X//////6H ERROR//31H VECTORS DEFINING THE PLANE ARE,      
     & 17H LINEAR DEPENDENT/34H ALL POINTS LIE ON A STRAIGHT LINE       
     & //13H TAKE ANOTHER,14H SET OF POINTS)                            
   12 FORMAT(11H0   A=     ,3F10.4,10H   B=     ,3F10.4)                
   13 FORMAT(28H0 VECTORS SPANNING THE PLANE)                           
   14 FORMAT(16H0 CELL CONSTANTS//5X,1HA,9X,1HB,9X,1HC,26X,             
     &5HALPHA,6X,4HBETA,7X,5HGAMMA//1X,3F10.4,2X,10HANGSTROEMS,5X,      
     &3F11.4,2X,7HDEGREES)                                              
   15 FORMAT(21H0 ORIGIN OF THE PLANE//1X,3F10.4)                       
   16 FORMAT(34H0 AREA OF THE PLANE LENGTH * WIDTH//1X,F10.4            
     &,5H    *,F10.4)                                                   
   17 FORMAT(1H0,14H METRIC MATRIX//3(1X,3F10.4//),8H  POINTS           
     &,19H DEFINING THE PLANE//30H     IN FRACTIONAL COORDINATES,19X,   
     &19HIN GRID COORDINATES//1X,3(3F10.4,15X,2F10.4,7X,A8,2X,6HPOINT   
     & //1X))                                                           
   18 FORMAT(23H0 PLOT OF DEFINED PLANE//22(3X,80A1/))                  
   19 FORMAT(28H1PARAMETERS OF SPECIAL PLANE/1X,27(1H-)//)              
   21 FORMAT(1H0)                                                       
C                                                                       
      WRITE(IOUT,19)                                                    
C                                                                       
      IF(IERR.EQ.0)GOTO 290                                             
      WRITE(IOUT,11)                                                    
      CALL EXIT                                                         
C                                                                       
  290 PP1(1)=0.                                                         
      PP1(2)=0.                                                         
      PP2(1)=FLOAT(IL)                                                  
      PP2(2)=0.                                                         
      PP3(2)=FLOAT(ILY)                                                 
      IF(ITYPE.EQ.1)PP3(2)=FLOAT(ILY)*0.5                               
      PP3(1)=C*FLOAT(IL)                                                
C                                                                       
      ALNU(1,1)=AI                                                      
      DO 300 N=1,21                                                     
      DO 300 M=2,80                                                     
  300 ALNU(M,N)=AB                                                      
      DO 400 N=1,22                                                     
  400 ALNU(1,N)=AI                                                      
      DO 500 M=1,80                                                     
  500 ALNU(M,22)=AMI                                                    
      ALNU(3,1)=AYY                                                     
      ALNU(5,1)=AA                                                      
      ALNU(6,1)=AXX                                                     
      ALNU(7,1)=AI                                                      
      ALNU(8,1)=ASS                                                     
      ALNU(71,21)=AXX                                                   
      ALNU(73,21)=AA                                                    
      ALNU(74,21)=AXX                                                   
      ALNU(75,21)=AI                                                    
      ALNU(76,21)=ASS                                                   
      II=3                                                              
      IF(ITYPE.EQ.1)II=11                                               
      IF(C.LT.1)GOTO 600                                                
      DO 510 N=3,19                                                     
      IP2=IFIX(60./C)+10                                                
      ALNU(IP2,N)=AS                                                    
  510 ALNU(10,N)=AS                                                     
      DO 530 N=10,IP2                                                   
      ALNU(N,3)=AS                                                      
  530 ALNU(N,19)=AS                                                     
      ALNU(70,II)=AP                                                    
      ALNU(71,II)=AP3                                                   
      ALNU(10,19)=AP                                                    
      ALNU(11,19)=AP1                                                   
      ALNU(12,19)=AB                                                    
      ALNU(IP2,19)=AP                                                   
      ALNU(IP2+1,19)=AP2                                                
      ALNU(IP2-1,19)=AB                                                 
      GOTO 800                                                          
  600 IF(CD.GT.0.)GOTO 700                                              
      IP1=IFIX(60./(1.-C))                                              
      IP1=70-IP1                                                        
      ALNU(10,II)=AP                                                    
      ALNU(11,II)=AP3                                                   
      DO 610 N=IP1,70                                                   
      ALNU(N,19)=AS                                                     
  610 ALNU(N,3)=AS                                                      
      DO 620 N=3,19                                                     
      ALNU(IP1,N)=AS                                                    
  620 ALNU(70,N)=AS                                                     
      ALNU(IP1,19)=AP                                                   
      ALNU(IP1+1,19)=AP1                                                
      ALNU(IP1+2,19)=AB                                                 
      ALNU(70,19)=AP                                                    
      ALNU(71,19)=AP2                                                   
      ALNU(69,19)=AB                                                    
      GOTO 800                                                          
  700 IP3=IFIX(C*60.)+10                                                
      DO 710 N=10,70                                                    
      ALNU(N,3)=AS                                                      
  710 ALNU(N,19)=AS                                                     
      DO 720 N=3,19                                                     
      ALNU(10,N)=AS                                                     
  720 ALNU(70,N)=AS                                                     
      ALNU(IP3-1,II)=AB                                                 
      ALNU(IP3,II)=AP                                                   
      ALNU(IP3+1,II)=AP3                                                
      ALNU(IP3+2,II)=AB                                                 
      ALNU(10,19)=AP                                                    
      ALNU(11,19)=AP1                                                   
      ALNU(12,19)=AB                                                    
      ALNU(70,19)=AP                                                    
      ALNU(71,19)=AP2                                                   
      ALNU(69,19)=AB                                                    
C                                                                       
      IF(IA.EQ.0)GOTO800                                                
      XX=PP(1)/FLOAT(IX)*60.+10.                                        
      IP1=IFIX(XX)                                                      
      YY=PP(2)/FLOAT(IY)*19.+3.                                         
      IP2=IFIX(YY)                                                      
      ALNU(IP1-2,IP2)=AA                                                
      ALNU(IP1-1,IP2)=ATT                                               
      ALNU(IP1,IP2)=AOO                                                 
      ALNU(IP1+1,IP2)=AMM                                               
C                                                                       
  800 WRITE(IOUT,14)AL,BL,CL,CAL,CBE,CGA                                
      WRITE(IOUT,17)G,P1,PP1,NAME(1),P2,PP2,NAME(1),P3,PP3,             
     & NAME(ITYPE+1)                                                    
      WRITE(IOUT,15)X,Y,Z                                               
      WRITE(IOUT,16)ALANG,BLANG                                         
      WRITE(IOUT,13)                                                    
      WRITE(IOUT,12)A,AM                                                
      IF(MP.LT.70)WRITE(IOUT,21)                                        
      WRITE(IOUT,18)ALNU                                                
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE CORD(A,B,IXI,IXF,IYI,IYF,I,X,Y,Z,PIJ1,PIJ2,PIJ3)       
C                                                                       
C     SUBROUTINE CALCULATES COORDINATES OF SPECIAL PLANE  DEFINED       
C                       BY ROUTINE PLANE.                               
C                                                                       
C              BY ROLF BACHMANN    DECEMBER 1980                        
C                                                                       
C                                                                       
C     INPUT   UNIT VECTORS A, B; THE GRID (IXF-IXI)*IYF-IYI)            
C             THE  ORIGIN  OF  THE  PLANE  X,Y,Z; AN INDEX I            
C     OUTPUT  THE COORDINATES OF  POINT PIJ PIJ1, PIJ2, PIJ3            
C                                                                       
      DIMENSION A(3), B(3)                                              
C                                                                       
C---- NUMBER OF POINTS IN X, Y DIRCTION                                 
C                                                                       
      IX=IXF-IXI                                                        
      IY=IYF-IYI                                                        
C                                                                       
C                                                                       
C---- CALCULATES GRID COORDINATES                                       
C                                                                       
      IX1=MOD(I,(IX+1))                                                 
      C=FLOAT(I)/FLOAT(IX+1)                                            
      IY1=INT(C)                                                        
      AX=FLOAT(IX1)                                                     
      AY=FLOAT(IY1)                                                     
C                                                                       
C---- THE COORDINATES IN FRACTIONS                                      
C                                                                       
      PIJ1=X+AX*A(1)+AY*B(1)                                            
      PIJ2=Y+AX*A(2)+AY*B(2)                                            
      PIJ3=Z+AX*A(3)+AY*B(3)                                            
C                                                                       
C---- INDEX=INDEX+1                                                     
C                                                                       
      I=I+1                                                             
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      FUNCTION XGY(X,G,Y)                                               
C                                                                       
      DIMENSION X(3),Y(3),G(3,3)                                        
C                                                                       
C ----SCALAR PRODUCT OR TWO VECTORS                                     
C                                                                       
      XGY=X(1)*(G(1,1)*Y(1)+G(1,2)*Y(2)+G(1,3)*Y(3))+                   
     &    X(2)*(G(2,1)*Y(1)+G(2,2)*Y(2)+G(2,3)*Y(3))+                   
     &    X(3)*(G(3,1)*Y(1)+G(3,2)*Y(2)+G(3,3)*Y(3))                    
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ATOPLA(P1,P2,P3,ATO,IX,IY,G,PP)                        
C                                                                       
C     PROGRAMME FOR  DETERMINIG THE  POSITION OF AN ATOM WITHIN         
C     A SPECIAL PLANE                                                   
C                                                                       
C                BY ROLF BACHMANN    DECEMBER 1980                      
C                                                                       
C     INPUT    THREE  POINTS  DEFINING  THE PLANE  P1,P2,P3,THE         
C              FRACTIONAL  COORDINATES  OF THE ATOM,  THE  GRID         
C              IX*IY AND THE METRIC G                                   
C     OUTPUT   A VECTOR PP GIVING  THE GRID  COORDINATES OF THE         
C              ATOM AND ITS DISTANCE  FROM THE PLANE                    
C                                                                       
      DIMENSION G(3,3),P1(3),P2(3),P3(3),ATO(3),E(3),C(3),PP(3)         
      DIMENSION A(3)                                                    
C                                                                       
      COMMON /APLANE/B(3)                                               
C                                                                       
C ----VECTOR C=ATO-P1,A=P2-P1                                           
C                                                                       
      A(1)=P2(1)-P1(1)                                                  
      A(2)=P2(2)-P1(2)                                                  
      A(3)=P2(3)-P1(3)                                                  
      C(1)=ATO(1)-P1(1)                                                 
      C(2)=ATO(2)-P1(2)                                                 
      C(3)=ATO(3)-P1(3)                                                 
C                                                                       
C ----NORMAL VECTOR OF THE PLANE                                        
C                                                                       
      PP(1)=XGY(A,G,C)/(XGY(A,G,A))                                     
      PP(2)=XGY(B,G,C)/(XGY(B,G,B))                                     
      E(1)=C(1)-PP(1)*A(1)-PP(2)*B(2)                                   
      E(2)=C(2)-PP(1)*A(2)-PP(2)*B(2)                                   
      E(3)=C(3)-PP(1)*A(3)-PP(2)*B(3)                                   
C                                                                       
C ----GRID COORDINATES OF THE ATOM AND DISTANCE TO THE PLANE            
C                                                                       
      PP(1)=XGY(A,G,C)/XGY(A,G,A)*FLOAT(IX)                             
      PP(2)=XGY(B,G,C)/XGY(B,G,B)*FLOAT(IY)                             
      PP(3)=SQRT(XGY(E,G,E))                                            
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE RANVEC(NPAR,RANR,seed,twopi)
C---- SUBROUTINE TO PROVIDE RANDOM VECTOR
C     NPAR number of elements
C     RANR random vector
C
      DIMENSION RANR(NPAR+1)
      INTEGER*4 SEED
C
C---- MAKE SURE TO GET ENOUGH RANDOM NUMBERS
C
      NPA = NPAR + 1
      DO 100 I=1,NPA
      RANR(I) = RAN(SEED)
  100 CONTINUE
c      type 99, ranr
   99 format (' ranr=',8f8.4)
C
C---- TRANSFORM TO STANDARD NORMAL RANR
C
      DO 200 I=1,NPAR,2
      RANSQ= SQRT(-2.*ALOG(RANR(I)))
      RANR(I) =  RANSQ*COS(TWOPI*RANR(I+1))
      RANR(I+1) = RANSQ*SIN(TWOPI*RANR(I+1))
  200 CONTINUE
C     type 99, ranr
      RETURN
      END
      SUBROUTINE RANCOR(NPAR,nij,cl,RANR,ranc)
C
C---- CALCULATE CORRELATED RANDOM VECTOR RANC
C
C     NPAR length of random vector
C
      DIMENSION cl(nij),RANR(NPAR+1),RANC(npar+1)
C
      do 10 i=1,npar+1
   10 ranc(i) = 0.0
      ij = 1
      do 200 j=1,npar
      do 100 i=1,j
      ranc(j) = ranc(j) + cl(ij)*ranr(i)
      ij = ij + 1
  100 continue
  200 continue
C      type 98, (cl(i),i=1,nij)
C   98 format(' cl = ',8f8.4)
C      type 99, ranc
C   99 format(' ranc= ',8f8.4)
      RETURN
      END
      subroutine crack(jirmin,jirmax,mcfp,vcv,c,cl,mdegn,mrank)
c
c---- crack var-covar for atom
c
      common  /vaco/ cv(40000),bffd(200)
      dimension vcv(mcfp,mcfp),c(mcfp,mcfp)
      dimension cl(3570)
      idif = mrank - jirmax
      ij = mdegn + 1
      do 10 i=1,mrank
      if(i.eq.jirmin) goto 11
      ij = ij + mrank - i + 1
   10 continue
   11 continue
      ii = 1
c     type 99, mdegn,mrank,jirmin,jirmax
c  99 format(' mdegn,mrank,jirmin,jirmax = ',4i5)
      do 200 i=jirmin,jirmax
      jj = ii
      do 100 j = i,jirmax
      vcv(jj,ii) = cv(ij)
C      vcv(jj,ii) = cv(ij)/(bffd(i)*bffd(j))
C      vcv(jj,ii) = sqrt(abs(vcv(jj,ii)))*sign(1.,vcv(jj,ii))
      jj = jj + 1
  100 ij = ij + 1
      ii = ii + 1
  200 ij = ij + idif
      call cormat(mcfp,vcv,c)
c
c---- transform c into linear array
c
      ij = 1
      do 400 i=1,mcfp
      do 300 j=1,i
      cl(ij) = c(i,j)
  300 ij = ij + 1
  400 continue
      return
      end
      SUBROUTINE CORMAT(NPAR,vcv,C)
C
C---- CALCULATE LOWER TRIANGULAR MATRIX C BY THE SQRT METHOD
C
C     NPAR number of elements of random vector(=rows in matrix)
C     vcv variance-covariance matrix of atom
C
      DIMENSION vcv(NPAR,NPAR),C(NPAR,NPAR)
C
C---- INITIALIZE
C
      DO 50 I=1,NPAR
      DO 50 J=1,NPAR
      C(I,J) = 0.0
   50 CONTINUE
c
c---- fill var-covar matrix
      do 90 i=1,npar
      do 90 j=i,npar
   90 vcv(i,j) = vcv(j,i)
C      type 98, ((vcv(j,i),j=i,npar),i=1,npar)
C
C---- CALCULATE MATRIX C
C
      SQ11 = SQRT(VCV(1,1))
      do 1000 i=1,npar     
      C(I,1) = VCV(I,1)/SQ11
      if(i.le.1) goto 1000
      if(i.eq.2) goto 3000
      do 2000 j=2,i-1
      ciksum = 0.
      do 300 k=1,i-1
      ciksum = ciksum + c(i,k)*c(j,k)
  300 continue
      c(i,j) = (vcv(i,j)-ciksum)/c(j,j)
 2000 continue
 3000 ciksum = 0.
      do 200 k=1,i-1
      ciksum = ciksum + c(i,k)**2
  200 continue
      c(i,i) = sqrt(abs(vcv(i,i)-ciksum))
 1000 continue
C      type 98, ((c(j,i),j=i,npar),i=1,npar)
C   98 format(' c(i,j)= ',8f9.5)
      RETURN
      END

      
  	      			
      

      
  	      			
      subroutine modify(mcfp,jirmin,jirmax,iatom,rc)
      dimension rc(mcfp+1),cv(40000),jir(84,60),bffd(200)
      dimension isot(60),occa(60),xyz(3,60),beta(6,60),
     &cr(10,60),dr(15,60),er(21,6),fr(28,6),ipara(1000),
     &occao(60),xyzo(3,60),betao(6,60),cro(10,60),dro(15,60),
     &ero(21,6),fro(28,6)
      common   /vaco/cv,bffd
      common   /parms/isot,occa,xyz,beta,cr,dr,er,fr,jir,ipara
      common   /parmol/occao,xyzo,betao,cro,dro,ero,fro
c
      j = 0
c     type 99, iatom
c  99 format (' iatom = ',i8)
      do 150 i=1,84
      if((jir(i,iatom)).eq.0) goto 150
      j = j + 1
      ip = jir(i,iatom)
c      bffd(i) = 1.0
      if(i.lt.2) goto 112
      if(i.lt.5) goto 116
      if(i.lt.11) goto 120
      if(i.lt.21) goto 124
      if(i.lt.36) goto 128
      if(i-57) 1281,1282,1282
  112 occa(iatom) = occao(iatom) + rc(j)
      goto 150
  116 xyz(i-1,iatom) = xyzo(i-1,iatom) + rc(j)
      goto 150
  120 beta(i-4,iatom) = betao(i-4,iatom) + rc(j)
      goto 150
  124 cr(i-10,iatom) = cro(i-10,iatom) + rc(j)
      goto 150
  128 dr(i-20,iatom) = dro(i-20,iatom) + rc(j)
      goto 150
 1281 er(i-35,iatom) = ero(i-35,iatom) + rc(j)
      goto 150
 1282 fr(i-56,iatom) = fro(i-56,iatom) + rc(j)
  150 continue
      return
      end

                       SUBROUTINE RESET(NA,KPARA,ndegn)                        
C                                                                       
C     THIS SUBROUTINE IS  ENTERED (FROM LOADNW)  ONCE FOR  SPECIAL PARA-
C     METERS AND ONCE FOR EACH ATOM TO COMPUTE THE UPDATED VALUES OL ALL
C     DEPENDENT PARAMETERS. LINEAR RELATIONSHIPS AND  RIGID-BODY PARAME-
C     TERS ARE COMPUTED DIRECTLY.    OTHER NON LINEAR  SYSYTEMS OF  CON-
C     STRAINTS MAY BE HANDELED BY THE USE THROUGH  SUBROUTINE * SPVAL *,
C                       SUPPLIED AS A DUMMY.                            
C                                                                       
C                     WRITTEN  BY L. W. FINGER                          
C                                                                       
C       ANHARMONIC THERMAL PARAMETERS OF FIFTH AND SIXTH ORDER ADDED.   
C                                                                       
C             FEBRUARY 1980  BY  U. ZUCKER  MPI-FKF STUTTGART.          
C                                                                       
      integer kdegn(1000)
      DIMENSION SAVE(60),ipara(1000)
C                                                                       
      common /con/ cnstr(1000)
      COMMON/parms/isot(60),occa(60),xyz(3,60),BETA(6,60),CR(10,60),
     &DR(15,60),er(21,6),fr(28,6),jir(84,60),ip(1000)
      equivalence (cnstr,kdegn)
C                                                                       
C
      I1(K) = MOD(K/256**3,256)                                         
      I2(K) = MOD(K/256**2,256)                                         
      I3(K) = MOD(K/256,256)                                            
      I4(K) = MOD(K,256)                                                
C                                                                       
C                                                                       
      do 3 i=1,1000
    3 ipara(i) = ip(i)
      LDEGN=NDEGN
C      type 98,kpara,ldegn
C   98 format(' kpara: ',i4,'   ldegn: ',i4)
C      type 97,(kdegn(j),cnstr(j+1),j=1,ldegn,2)
C   97 format(' cnstr: ',5(i10,f8.5))                                                     
C      type 96,(ipara(kpara+j),j=1,20)
C   96 format(' ipara: ',8i8)
C                                                                       
      IF(NA.GT.0)GO TO 100                                              
c     if(na.eq.0) type 9200
 9200 format(' NO SPECIAL PARAMATERS ALLOWED. CALL EXIT !')
      call exit
C---------------------------------------------------------------------- 
C                                                                       
C---- COMPUTE NEW VALUES OF PARAMETERS FOR ATOM NA.                     
C                                                                       
 100  NPARA=5                                                           
      IF(ISOT(NA).GE.1)NPARA=NPARA+5                                    
      IF(ISOT(NA).GE.2)NPARA=NPARA+10                                   
      IF(ISOT(NA).GE.3)NPARA=NPARA+15                                   
      IF(ISOT(NA).GE.4)NPARA=NPARA+21                                   
      IF(ISOT(NA).GE.5)NPARA=NPARA+28                                   
      IF(LDEGN.LE.0)GO TO 1005                                          
C                                                                       
      DO 1000 ID = 1,NPARA                                              
      J = ID                                                            
C                                                                       
C---- IS PARAMETER INDEPENDENT?                                         
C                                                                       
      IF(IPARA(KPARA+J).GT.0)GO TO 1000                                 
C                                                                       
C---- IS THIS PARAMETER DEPENDENT ON ANY OTHERS?                        
C                                                                       
      DO 900 K=1,LDEGN,2                                                
      IF(I3(kdegn(K)).NE.NA.OR.I4(kdegn(K)).NE.J)GO TO 900              
C                                                                       
C---- DEPENDENT PARAMETER. SET IT.                                      
C                                                                       
      IF(IPARA(KPARA+J).LT.0)GO TO 900                                  
      IPARA(KPARA+J)=-1                                                 
      TEMPAR=0.                                                         
      KK=K                                                              

      IF(I1(kdegn(K)).NE.0.OR.I2(kdegn(K)).NE.0)GO TO 110               
      TEMPAR=cnstr(K+1)                                                    
      KK=KK+2                                                           
 110  IF(KK.GT.LDEGN)GO TO 170                                          
      DO 150 L=KK,LDEGN,2                                               
      IF(I3(kdegn(L)).NE.NA.OR.I4(kdegn(L)).NE.J)GO TO 150              
C                                                                       
C---- IS THE INDEPENDENT PARAMETER SPECIAL?                             
C                                                                       
      MK=I2(kdegn(L))                                                   
C      IF(I1(kdegn(L)).EQ.0) type 92,i3(kdegn(l)),i4(kdegn(l)),l,j,mk,mj
C   92 format(' i3 =',i4,' i4 =',i4,' l =',i4,' j =',i4,' mk =',i4,
C     &' mj =',i4)
c     if(i1(kdegn(l)).eq.0) goto 200
C                                                                       
C---- NO. LINEAR DEPENDENCIES ONLY. WHAT KIND OF PARAMETER?             
C                                                                       
c     MJ=I1(kdegn(L))                                                   
      mj = na
c     IF(MK.LT.2)GO TO 112
      if(mk.lt.2)goto 1000                                              
      IF(MK.LT.5)GO TO 116                                              
      IF(MK.LT.11)GO TO 120                                             
      IF(MK.LT.21)GO TO 124                                             
      IF(MK.LT.36)GO TO 128                                             
      IF(MK-57)1281,1282,1282                                           
 112  TEMPAR=TEMPAR+OCCA(MJ)*cnstr(L+1)                                    
      ITYPE=1                                                           
      GO TO 150                                                         
 116  TEMPAR=TEMPAR+XYZ(MK-1,MJ)*cnstr(L+1)                                
      ITYPE=2                                                           
      GO TO 150                                                         
 120  TEMPAR=TEMPAR+BETA(MK-4,MJ)*cnstr(L+1)                               
      ITYPE=3                                                           
      GO TO 150                                                         
 124  TEMPAR=TEMPAR+CR(MK-10,MJ)*cnstr(L+1)                                
      ITYPE=4                                                           
      GO TO 150                                                         
 128  TEMPAR=TEMPAR+DR(MK-20,MJ)*cnstr(L+1)                                
      ITYPE=5                                                           
      GOTO 150                                                          
 1281 TEMPAR=TEMPAR+ER(MK-35,MJ)*cnstr(L+1)                                
      ITYPE=6                                                           
      GOTO 150                                                          
 1282 TEMPAR=TEMPAR+FR(MK-56,MJ)*cnstr(L+1)                                
      ITYPE=7                                                           
 150  CONTINUE                                                          
 151  GO TO(152,154,156,158,160,162,164),ITYPE                          
 152  OCCA(NA)=TEMPAR                                                   
      GO TO 1000                                                        
 154  XYZ(J-1,NA)=TEMPAR                                                
      GO TO 1000                                                        
 156  BETA(J-4,NA)=TEMPAR                                               
      GO TO 1000                                                        
 158  CR(J-10,NA)=TEMPAR                                                
      GO TO 1000                                                        
 160  DR(J-20,NA)=TEMPAR                                                
      GO TO 1000                                                        
 162  ER(J-35,NA)=TEMPAR                                                
      GO TO 1000                                                        
 164  FR(J-56,NA)=TEMPAR                                                
      GO TO 1000                                                        
 170  IF(J.EQ.1) ITYPE=1                                                
      IF(J.GE.2) ITYPE=2                                                
      IF(J.GE.5) ITYPE=3                                                
      IF(J.GE.11)ITYPE=4                                                
      IF(J.GE.21)ITYPE=5                                                
      IF(J.GE.36)ITYPE=6                                                
      IF(J.GE.57)ITYPE=7                                                
      GO TO 151                                                         
C                                                                       
  200 type 9200
      call exit
c
c
  900 continue
 1000 CONTINUE                                                          
C1001 DO 1005 I=1,NPARA                                                 
C     IF(IPARA(KPARA+I).LT.0)IPARA(KPARA+I)=0                           
 1005 CONTINUE                                                          
C                                                                       
C      type 99,(cr(j,na),j=1,10)
C   99 format(' cr: ',10f8.5)
      RETURN                                                            
      END                                                               
                       SUBROUTINE REPAK(natom,mdegn)                        
C                                                                       
C                                                                       
      integer kdegn(1000)
      dimension ipara(1000)
C                                                                       
      common /con/ cnstr(1000)
      COMMON/parms/isot(60),occa(60),xyz(3,60),BETA(6,60),CR(10,60),
     &DR(15,60),er(21,6),fr(28,6),jir(84,60),ip(1000)
      equivalence (cnstr,kdegn)
C                                                                       
C
      I1(K) = MOD(K/256**3,256)                                         
      I2(K) = MOD(K/256**2,256)                                         
      I3(K) = MOD(K/256,256)                                            
      I4(K) = MOD(K,256)                                                
C                                                                       
C               
      do 100 i=1,mdegn,2
      ind1 = i1(kdegn(i))
      ind2 = i2(kdegn(i))
      id1 = i3(kdegn(i))
      id2 = i4(kdegn(i))
      do 200 j=1,natom
      do 200 k=1,84
      if(jir(k,j).eq.ind2) goto 250
  200 continue
  250 kdegn(i) = ipak(j,k,id1,id2)
c     type 99,j,k,id1,id2
c  99 format(' i1 =',i4,' i2 =',i4,' i3 =',i4,' i4 =',i4)
  100 continue
      return
      end
C
      function ipak(i1,i2,i3,i4)
      ipak = ((i1*256 + i2) * 256 + i3) * 256 + i4
      return
      end
                                                        
      subroutine mocad(izaehl,isq)
      dimension qp(10000),qm(10000)
      dimension array(10000),qs(10000),iplus(10000),iminus(10000)
      common /bach/ array,qp,qm,qs,iplus,iminus 
      do 100 i=1,izaehl
      if (isq.eq.1) x = (array(i) - qs(i))**2
      if (isq.eq.0) x = qs(i)-array(i)
      if(x) 10,20,20
   10 qm(i) = qm(i) + x
      iminus(i) = iminus(i) + 1
      goto 100
   20 qp(i) = qp(i) + x
      iplus(i) = iplus(i) + 1
  100 continue
      return
      end
c
      subroutine mocen(izaehl,isq)
      dimension qp(10000),qm(10000)
      dimension array(10000),qs(10000),iplus(10000),iminus(10000)
      common /bach/ array,qp,qm,qs,iplus,iminus
      do 100 i=1,izaehl
      if(iminus(i).gt.0)qm(i) = qm(i)/float(iminus(i))
      if(iplus(i).gt.0)qp(i) = qp(i)/float(iplus(i))
      if (isq.eq.1) qp(i) = sqrt(qp(i))
c      qm(i) = -sqrt(abs(qm(i)))
C      qm(i) = qs(i) - qm(i)
C      qp(i) = qs(i) - qp(i)
  100 continue
      return
      end
     
