C************************************************************* C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C************************************************************** C PROGRAM HKLPLOT implicit none INTEGER ISIZ,NPARM,MCOLS PARAMETER (ISIZ= 2346500) PARAMETER (NPARM=200) PARAMETER (MCOLS=500) EXTERNAL SETDEV INTEGER SETDEV C Hint to linker for block data: EXTERNAL HKBLK C C---- HKLPLOT file variables C DIMENSION RSYM(4,4,192) CHARACTER*80 TITLE CHARACTER*10 VERSN,PGNAME,NAMSPG,LTYPE*1 LOGICAL SQUARE,CALFLG,HPFLG,LYRESO REAL RSYM,RLIM2,HKLPDB,FMXX,SIG,ANOM,FSCAL,RR,RO,RF REAL RESMIN,RESMAX,SMIN,SMAX,RLIM,PROJ,PRJLIM,ARRAY(ISIZ) INTEGER I,ITOK,NSYM,NSYMB,NTERM,NSPGRP,NMULT,IHPROJ(3) INTEGER LIST,IGEN,NLPRGI,LDUM,IFAIL,NSSQ4,NSMP,NERR,IFP INTEGER NL,NREFL,NSYMP,NCOLS,NREFLN,IHMAX,IKMAX,ILMAX, + NCODE,J,IHMX21,IKMX21,ICHK,ICHK2 REAL RLIM12,RLIM22,RLIM32,VOL C C---- MTZ file variables C INTEGER LOOKUP,MTZIN REAL CELLMTZ,CELL,RANGES,RCELL,RVOL DIMENSION LOOKUP(MCOLS),CELLMTZ(6),CELL(6) DIMENSION RANGES(2,MCOLS),RCELL(6) C C---- things for parser C CHARACTER KEY*4,CVALUE(NPARM)*4,LINE*80,HKLZON(20)*3 REAL FVALUE(NPARM) LOGICAL LEND,LPRINT INTEGER LUNIN,LUNOUT,IPRINT INTEGER NTOK,IBEG(NPARM),IEND(NPARM),ITYPE(NPARM),IDEC(NPARM) CHARACTER CTPRGI(MCOLS)*1,LSPRGI(MCOLS)*30 C C---- Common Blocks C COMMON /SYMMM/ NSYM,NMULT,RSYM COMMON /PINOUT/ LUNIN,LUNOUT COMMON /HKLPLTC/ TITLE,HKLZON COMMON /HKLPLT1/ IFP,SQUARE,CALFLG,HPFLG C Note - there is no way to set NTERM - always = 0 now. C Redundant code should be removed - some day. COMMON /HKLPLT2/ NSYMB,SMIN,SMAX,RLIM,RLIM2,NTERM,NSPGRP, + HKLPDB, FMXX,SIG,IHPROJ,PROJ,PRJLIM,LIST,ANOM,FSCAL,IGEN COMMON /HKLPLT3/RR(3,3,6) C C---- cell orthogonalisation variables C COMMON /ORTHOG/ RO(4,4),RF(4,4) C C---- although this common really belongs to rwbrook it is useful to C keep it here. C C C---- Data Statements C DATA NLPRGI,LSPRGI/6,'H','K','L','F','SIGF','DANO', + 494*' '/ DATA CTPRGI /MCOLS*' '/ C DATA CTPRGI/'H','H','H','F','Q','D', C + 494*' '/ C LUNIN=5 LUNOUT=6 CALL CCPFYP CALL MTZINI CALL CCPDPN(LUNIN,'DATA','READONLY','F',LDUM,IFAIL) CALL CCPDPN(LUNOUT,'PRINTER','UNKNOWN','F',LDUM,IFAIL) CALL CCPRCS (LUNOUT,'HKLPLOT','$Date: 2002/08/05 14:06:47 $') HKLZON(1) = 'FIN' CELL(1) = 0.0 FSCAL = 1.0 HKLPDB = 0 NTERM = 0 IGEN = 0 ANOM = 0 SIG =0 LIST =0 NSSQ4=100 IHPROJ(1) = 0 IHPROJ(2) = 0 IHPROJ(3) = 0 PROJ = 0 RLIM = 2.0 SMIN = 0.00 SMAX = 0.25 LYRESO = .FALSE. NSYMB=6 NSMP=0 NERR=0 NMULT =1 1 CONTINUE NTOK=NPARM LINE=' ' CALL PARSER (KEY,LINE,IBEG,IEND,ITYPE,FVALUE,CVALUE,IDEC,NTOK, + LEND,.TRUE.) IF (LEND) GO TO 100 IF (KEY.EQ.'TITL') THEN TITLE=LINE(IBEG(2):80) C C---- HKLPDB C ELSE IF (KEY.EQ.'HKLP') THEN WRITE (LUNOUT,'(2x,a)') + ' HKLPDB style file will be output for display in FRODO or O' HKLPDB = 1 C C C---- HKLZONes C ELSE IF (KEY.EQ.'ZONE') THEN IF(NTOK.GT.20) NTOK = 20 DO 2 ITOK = 2,NTOK HKLZON(ITOK-1) = LINE(IBEG(ITOK):IBEG(ITOK)+2) 2 CONTINUE HKLZON(NTOK) = 'FIN' WRITE (LUNOUT,'(20(2x,a))') ' ZONES',(HKLZON(I),I=1,NTOK) C C C---- SYMBOL number _ see GINO writeup (diamond = default) C ELSE IF (KEY.EQ.'SYMB') THEN CALL GTPINT(2,NSYMB,NTOK,ITYPE,FVALUE) WRITE (LUNOUT,8001) NSYMB 8001 FORMAT (///,' SYMBOL no - see GINO writeup diamond default',I5) C C---- resolution limits - either 4SINSQ/LANDASQ or D C ELSE IF (KEY.EQ.'RESO') THEN ITOK = 2 C C ***************************************************** CALL RDRESO(ITOK,ITYPE,FVALUE,NTOK,RESMIN,RESMAX,SMIN,SMAX) C ***************************************************** LYRESO = .TRUE. C c ELSE IF (KEY.EQ.'SQUA') THEN 40 WRITE (LUNOUT,*) ' Input FS will be squared' SQUARE= .TRUE. c c ELSE IF (KEY(1:3).EQ.'BIN') THEN CALL GTPINT(2,NSSQ4,NTOK,ITYPE,FVALUE) IF(NSSQ4 .GT.1000) NSSQ4 = 1000 WRITE (LUNOUT,*) ' Resetting number of bins to:',NSSQ4 c ELSE IF (KEY(1:3).EQ.'DEV') THEN KEY = LINE(IBEG(2):IEND(2)) NTERM = SETDEV (KEY) IF (NTERM.LT.0) CALL CCPERR(1, ' Device not set ') c C---- SYMM or NSPGRP ELSE IF (KEY.EQ.'NSPG'.OR.KEY.EQ.'SYMM') THEN ITOK = 2 C C *************************************************** CALL RDSYMM(ITOK,LINE,IBEG,IEND,ITYPE,FVALUE,NTOK,NAMSPG,NSPGRP, + PGNAME,NSYM,NSMP,RSYM) IPRINT = 1 LPRINT = .TRUE. CALL PGDEFN(PGNAME,NSMP,NSYM,RSYM,LPRINT) C *************************************************** C NMULT = NSYM/NSMP C---- we need to set up centric and epsilon tests C CALL EPSLN (NSYM,NSMP,RSYM,IPRINT) C NSYM = NSMP C WRITE (LUNOUT,8003)NSPGRP 8003 FORMAT (/,' Spacegroup No',I6) WRITE (LUNOUT,8004)NMULT 8004 FORMAT (/,' Lattice multiplicity' ,I6) C C---- sigma cutoff for reflections C ELSE IF (KEY(1:3).EQ.'SIG') THEN CALL GTPREA(2,SIG,NTOK,ITYPE,FVALUE) WRITE (LUNOUT,8005)SIG 8005 FORMAT (///,' Sigma cutoff for reflections ',F8.3) C C---- h k l etc for projection C ELSE IF (KEY.EQ.'PROJ') THEN CALL GTPINT(2,IHPROJ(1),NTOK,ITYPE,FVALUE) CALL GTPINT(3,IHPROJ(2),NTOK,ITYPE,FVALUE) CALL GTPINT(4,IHPROJ(3),NTOK,ITYPE,FVALUE) CALL GTPREA(5,PROJ,NTOK,ITYPE,FVALUE) CALL GTPREA(6,PRJLIM,NTOK,ITYPE,FVALUE) IF (PRJLIM.EQ.0.0) PRJLIM = 0.002 WRITE (LUNOUT,9005)IHPROJ,PRJLIM,PROJ 9005 FORMAT (///, + ' Intensities plotted perp to this axis',3I4,/, + ' - all reflections LT ',F8.4,' from reciprocal height',F8.4, + ' plotted') C C---- anom plot - plot pos or neg ie F(+) or F(-) C ELSE IF (KEY.EQ.'ANOM') THEN ANOM = 0.5 KEY= LINE(IBEG(2):IBEG(2)) IF (KEY(1:1).EQ.'P' .OR. KEY(1:1).EQ.'p') ANOM = 0.5 IF (KEY(1:1).EQ.'N' .OR. KEY(1:1).EQ.'n') ANOM = -0.5 C IF (ANOM.EQ.0.5) WRITE (LUNOUT,'(//,A)') + ' Anomalous F(+) value plotted - assign DANO for HKLIN' C IF (ANOM.EQ.-0.5) WRITE (LUNOUT,'(//,A,/,A)') + ' Anomalous F(-) value plotted - assign DANO for HKLIN', + ' the axes are still labelled H K L >= 0 but the plotted', + ' F = F(-)' C WRITE (LUNOUT,'(//,A,//)') + ' NB - you should have a non-centric spacegroup number!' C C---- list H K L which are plotted C ELSE IF (KEY.EQ.'LIST') THEN LIST = 1 WRITE (LUNOUT,'(/,A,/)') ' Plotted reflections listed' C C---- fscale scale applied to SYMBOLs on output C ELSE IF (KEY.EQ.'SCAL') THEN IF(ITYPE(2).EQ.2)CALL GTPREA(2,FSCAL,NTOK,ITYPE,FVALUE) IF(ITYPE(3).EQ.2)CALL GTPREA(3,FSCAL,NTOK,ITYPE,FVALUE) WRITE (LUNOUT,'(/,A,F6.2/)') ' Output SYMBOLs scaled by',FSCAL C C---- cell - you can use this to generate a reciprocal lattice pattern C for a given cell C ELSE IF (KEY.EQ.'CELL') THEN IGEN = 1 CALL RDCELL(2,ITYPE,FVALUE,NTOK,CELL) WRITE (LUNOUT,'(/,A,/,10X,6F8.2/)') + ' Input cell used for reciprocal lattice generation ',CELL C---- "LABI" i.e. set input labels C ELSE IF (KEY.EQ.'LABI') THEN C C ***************************************** CALL LKYIN(1,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND) C ***************************************** C C Must have first 4 labels h k l F C C---- end of input check for errors C ELSE IF (KEY(1:3).EQ.'END') THEN GO TO 100 ELSE WRITE (LUNOUT,'(/,A,3X,A,/)')' Unrecognized keyword :',KEY NERR=NERR+1 ENDIF GO TO 1 C C input finished C 100 IF (NERR.GT.0) THEN WRITE (LUNOUT,9901) NERR 9901 FORMAT (/1X,I6,' Errors in input, see above !!!!'/) Call ccperr(1,' Fatal error') ENDIF IF (NTERM.LT.0) Call ccperr(1, ' Device not set ') DO 201 I=1,4 LOOKUP(I)=-1 201 CONTINUE C If cell given - assume you want to generate a pattern IF(CELL(1).GT.0.0) GO TO 210 C No need to read S - workit out if not there. LOOKUP(5) = 0 LOOKUP(6)=0 IGEN = 0 NL=100 NREFL=0 C C---- this read in the header information from an MTZ file MTZIN = 1 IPRINT = 1 C *********************************** CALL LROPEN(MTZIN,'HKLIN',IPRINT,IFAIL) C *********************************** C IF (IFAIL.EQ.1) CALL CCPERR(1,' No input file???????') C C ***************************************** CALL LRASSN(MTZIN,LSPRGI,NLPRGI,LOOKUP,CTPRGI) C ***************************************** C USE SIG as a SIGMA input flag - nonzero value means SIGF input IF(LOOKUP(5).NE.0 .AND. SIG.EQ. 0) SIG = -1.0 C C---- Then call lrcell etc... C C ********************* CALL LRCELL(MTZIN,CELLMTZ) C ********************* CELL(1) = CELLMTZ(1) CELL(2) = CELLMTZ(2) CELL(3) = CELLMTZ(3) CELL(4) = CELLMTZ(4) CELL(5) = CELLMTZ(5) CELL(6) = CELLMTZ(6) C IF (.NOT.LYRESO) THEN CALL LRRSOL(MTZIN,SMIN,SMAX) ENDIF C C C---- Read symmetry and check C C ********************************************** CALL LRSYMI(MTZIN,NSYMP,LTYPE,NSPGRP,NAMSPG,PGNAME) IF(NSMP.GT.0 .AND.NSYMP.NE.NSMP) CALL CCPERR + (1,' Disagreement between symmetry operators') CALL LRSYMM(MTZIN,NSYM,RSYM) C ********************************************** C Find largest FP CALL LRINFO(MTZIN,VERSN,NCOLS,NREFLN,RANGES) FMXX = RANGES(2,LOOKUP(4)) C C NMULT = NSYM/NSYMP C---- we need to set up centric and epsilon tests C CALL EPSLN (NSYM,NSYMP,RSYM,IPRINT) C NSYM = NSYMP C WRITE (LUNOUT,8003)NSPGRP WRITE (LUNOUT,8004)NMULT C C C---- come here if cell given and no MTZ file input C 210 CONTINUE C SLIM = SMAX RLIM = 1./SQRT(SMAX) WRITE (LUNOUT,8002)SMAX,RLIM 8002 FORMAT (/,' 4SINSQ/LSQ ',F10.5,' resolution limit',F10.5) C WRITE (LUNOUT, '(/A,6F8.2)') ' Cell dimensions:',CELL IHMAX=INT(CELL(1)/RLIM) + 1 IKMAX=INT(CELL(2)/RLIM) + 1 ILMAX=INT(CELL(3)/RLIM) + 1 C C---- find rlim2 to allow for axis labels at position hmax+1 ,kmax+1 lmax+1 C RLIM12= CELL(1)/(IHMAX+1 -1) RLIM22= CELL(2)/(IKMAX+1 -1) RLIM32= CELL(3)/(ILMAX+1 -1) RLIM2=RLIM12 IF (RLIM22.LT.RLIM2)RLIM2=RLIM22 IF (RLIM32.LT.RLIM2)RLIM2=RLIM32 CALL RBFRO1(CELL,VOL,RR) C C---- set NCODE eq 1 to get an orthogonalising matrix to work out S C it will have to be reset later for different projections C NCODE=1 RO(4,4)=1 DO 44 I=1,3 RO(4,I)=0 RO(I,4)=0 DO 4 J=1,3 RO(J,I)=RR(J,I,NCODE) 4 CONTINUE 44 CONTINUE CALL RBRINV(RO,RF) IHMX21=2*IHMAX +1 IKMX21=2*IKMAX +1 ICHK= IHMX21*IKMX21*(ILMAX+1) ICHK2= 2*ICHK IF (2*ICHK.GT.ISIZ) THEN WRITE(LUNOUT,'(2A,I9,A,/,A,3I5)') + ' Array dimensions too small ', + ' - decrease resolution or increase ISIZ to ',ICHK2, + ' This equals ( 2*HMAX+1) * ( 2*KMAX+1) * (LMAX+1) *2', + ' Hmax Kmax Lmax are:',IHMAX,IKMAX,ILMAX Call ccperr + (1, ' Array dimensions too small - decrease resolution') END IF IF(HKLPDB .NE.0.0) THEN C Find reciprocal cell etc. CALL RBRCEL(RCELL,RVOL) CALL CCPDPN(20,'XYZOUT','NEW','F',0,IFAIL) C For pseudo PDB multiply reciprocal lattice by 1000.0 WRITE(20,'(A,3F9.3,3F7.2)' )'CRYSTL',RCELL c RCELL(1)=0.001*RO(2,1) c RCELL(2)=0.001*RO(2,2) c RCELL(3)=0.001*RO(2,3) c RCELL(4)=0.001*RO(2,4) C I need to write out RO transpose since HKL being output as C pseudo XO YO ZO ... WRITe(20,'(a,3F10.5,5X,3F10.5)') + 'SCALE1 ',(RO(J,1),J=1,4) WRITe(20,'(a,3F10.5,5X,3F10.5)') + 'SCALE2 ',(RO(J,2),J=1,4) WRITe(20,'(a,3F10.5,5X,3F10.5)') + 'SCALE3 ',(RO(J,3),J=1,4) END IF C 1000 CONTINUE CALL LCFPREC(ARRAY,IHMX21,IKMX21,ILMAX,NSSQ4) IF (HKLPDB .NE. 0) CLOSE (UNIT=20,STATUS='KEEP') CALL CCPERR(0,' End of job') C C---- Error Messages C 9990 Call ccperr(1, ' File definition error') 9000 Call ccperr(1, ' Input data must end with END') END C C SUBROUTINE BEGPLT(FILNAM) C ========================== C CHARACTER*(*) FILNAM INTEGER NTERM,NDEVIC,NDIREC C INTEGER LENSTR EXTERNAL LENSTR,GSINIT,GSDVIC,GSRFNT,GSPICT,GSPRNT C COMMON /DEVATT/ NTERM,NDEVIC,NDIREC C SAVE C C---- NTERM = 1 Terminal = 0 Not a Terminal C NDEVIC = 1 Undefined = 2 Paper = 3 VT640 C NDIREC = 0 Not Interactive = 1 Interactive C IF (LENSTR(FILNAM).EQ.0) THEN CALL GSINIT ('PLOT.PLT') ELSE CALL GSINIT (FILNAM) ENDIF CALL GSPRNT (0) IF (NTERM.EQ.1) CALL GSDVIC (NDIREC,NDEVIC) CALL GSRFNT CALL GSPICT END C C BLOCK DATA HKBLK CHARACTER TITLE*80, HKLZON(20)*3 COMMON /HKLPLTC/ TITLE,HKLZON SAVE DATA HKLZON/20*' '/ END C C SUBROUTINE CHRINT (INUM,IWIDTH) C ================================ C INTEGER INUM,IWIDTH,NJUST REAL SIZEX,SIZEY C COMMON /XYSIZ/ SIZEX,SIZEY C SAVE C DATA NJUST/2/ C CALL GSINUM (INUM,IWIDTH,SIZEX,SIZEY,NJUST) END C C SUBROUTINE CHRSIZ (XCHAR,YCHAR) C =============================== C REAL XCHAR,YCHAR,SIZEX,SIZEY C COMMON /XYSIZ/ SIZEX,SIZEY C SAVE C SIZEX = XCHAR SIZEY = YCHAR C CALL GSSCLC (XCHAR,YCHAR) END C C SUBROUTINE CHRSTR (STRING) C ========================== C CHARACTER*(*) STRING C CALL GSSTRC (STRING) END C C SUBROUTINE DASHLN (W,X,Y,Z) C =========================== C REAL W,X,Y,Z C CALL GSCOLR (3) END C C SUBROUTINE DEVICE (INUM) C ======================== C INTEGER INUM C IF (INUM.EQ.3) THEN CALL GSMIXC (1) ENDIF END C C SUBROUTINE DSYMB (INUM) C ======================= C INTEGER INUM,NSET REAL SIZEX,SIZEY C COMMON /XYSIZ/ SIZEX,SIZEY C SAVE C DATA NSET/1/ C CALL GSGSYS (INUM,NSET,0.0,0.0,SIZEX,SIZEY) END C C SUBROUTINE ENDPIC (IDUM) C ======================== C INTEGER IDUM C CALL GSENDP CALL GSPICT END C C SUBROUTINE ENDPLT (IDUM) C ======================== C INTEGER IDUM C CALL GSSTOP END SUBROUTINE GRAPHDRAW(ARRAY,IHMX21,IKMX21,ILMAX, + NREFL,HMAX,KMAX,LMAX,HMIN,KMIN,LMIN, + TOP,BOTTOM,LEFT,RIGHT,FMIN,FMAX,SECTION, + DIRECTION,RF,HKLMAX,DEVCEN,DEVSCL) C ============================================================= C C---- to plot on the screen the data supplied C REAL ARRAY(2,IHMX21,IKMX21,0:ILMAX) C C---- LCFPREC file variables C CHARACTER*2 AXSLAB(3) CHARACTER*80 TITLE,HKLZON(20)*3 LOGICAL SQUARE,CALFLG,HPFLG REAL SMIN,SMAX,RLIM,PROJ,PRJLIM INTEGER NSYMB,NTERM,NSPGRP,NMULT,INDLAB(3),IHPROJ(3) INTEGER NREFL,HMAX,KMAX,LMAX,HMIN,KMIN,LMIN, + HKLMAX,LUNIN,LUNOUT INTEGER SECTION,DIRECTION,J REAL TOP,BOTTOM,LEFT,RIGHT,EXTRA REAL POSX,POSY,FMIN,FMAX,DIF,SIZE,RF(4,4) REAL AHO,AKO,ALO REAL DEVCEN(0:4,2),DEVSCL(0:4) DIMENSION RSYM(4,4,192) C C---- Common Blocks C COMMON /PINOUT/ LUNIN,LUNOUT COMMON /SYMMM/ NSYM,NMULT,RSYM COMMON /HKLPLTC/ TITLE,HKLZON COMMON /HKLPLT1/ IFP,SQUARE,CALFLG,HPFLG COMMON /HKLPLT2/ NSYMB,SMIN,SMAX,RLIM,RLIM2,NTERM,NSPGRP, + HKLPDB, FMXX,SIG,IHPROJ,PROJ,PRJLIM,LIST,ANOM,FSCAL,IGEN SAVE /PINOUT/, /SYMMM/, /HKLPLTC/, /HKLPLT1/, /HKLPLT2/ C C---- Data Statement C DATA AXSLAB/'AS','BS','CS'/ C EXTRA=2.5 C C---- for SYMBOL size : C ejd try factor of 10 - still too small on hpa4 and calcomp C IF (CALFLG)EXTRA=EXTRA*10 IF (HPFLG)EXTRA=EXTRA*10 CALL CHRSIZ (1.0,1.0) DIF=FMAX-FMIN CALL MOVETOB (DEVCEN(NTERM,1),DEVCEN(NTERM,2)) CALL DSYMB (15) IMIN1 = HMIN IMAX1=HMAX INDLAB(1)=HMAX+1 IMIN2 = KMIN IMAX2=KMAX INDLAB(2)=KMAX+1 IMIN3=-LMAX IMAX3=LMAX INDLAB(3)=LMAX+1 IF (DIRECTION .EQ.1) THEN IMIN1=SECTION IMAX1=SECTION INDLAB(1)=SECTION ENDIF IF (DIRECTION .EQ.2) THEN IMIN2=SECTION IMAX2=SECTION INDLAB(2)=SECTION ENDIF IF (DIRECTION .EQ.3) THEN IMIN3=SECTION IMAX3=SECTION INDLAB(3)=SECTION ENDIF ISIGN=-1 DO 1 I3=IMIN3,IMAX3 IF (I3.GE.0)ISIGN=1 IND3=I3*ISIGN DO 21 I2=IMIN2,IMAX2 IND2=I2*ISIGN +KMAX +1 DO 2 I1=IMIN1,IMAX1 IND1=I1*ISIGN +HMAX +1 IF (ARRAY(1,IND1,IND2,IND3).LT.FMIN)GO TO 2 AHO=I1*RF(1,1)+I2*RF(2,1)+I3*RF(3,1) AKO=I1*RF(1,2)+I2*RF(2,2)+I3*RF(3,2) ALO=I1*RF(1,3)+I2*RF(2,3)+I3*RF(3,3) c CALL ORTMAT(RF,I1,I2,I3,AHO,AKO,ALO) C C---- projection axis - only interested in XO=PROJ C IF ((DIRECTION.EQ.-1) .AND. ABS(AHO-PROJ) .GT.PRJLIM) + GO TO 2 IF (LIST.EQ.1) THEN IFFF = NINT(ARRAY(1,IND1,IND2,IND3)) WRITE (LUNOUT,'(A,3I4,I12)') + ' Refln plotted', I1,I2,I3,IFFF ENDIF SIZE = FSCAL*(ARRAY(1,IND1,IND2,IND3)-FMIN)/DIF IF(SIZE.GT.1) SIZE = 1.0 IF (SQUARE)SIZE=SIZE*SIZE/2 6789 FORMAT (3I3,4F15.2,5X,3F15.2,I7) POSX=AKO*RLIM2*DEVSCL(NTERM)+DEVCEN(NTERM,1) POSY=ALO*RLIM2*DEVSCL(NTERM)+DEVCEN(NTERM,2) CALL MOVETOB (POSX,POSY) CALL CHRSIZ (SIZE*EXTRA,SIZE*EXTRA) CALL DSYMB (NSYMB) 60 CONTINUE 2 CONTINUE 21 CONTINUE 1 CONTINUE C C---- label a* b* c* C CALL CHRSIZ (2.5,2.5) I1=INDLAB(1) I2=0 I3=0 J=1 61 CONTINUE C CALL ORTMAT(RF,I1,I2,I3,AHO,AKO,ALO) AHO=I1*RF(1,1)+I2*RF(2,1)+I3*RF(3,1) AKO=I1*RF(1,2)+I2*RF(2,2)+I3*RF(3,2) ALO=I1*RF(1,3)+I2*RF(2,3)+I3*RF(3,3) POSX=AKO*RLIM2*DEVSCL(NTERM)+DEVCEN(NTERM,1) POSY=ALO*RLIM2*DEVSCL(NTERM)+DEVCEN(NTERM,2) CALL MOVETOB (POSX,POSY) IF (POSX.GT.LEFT.AND.POSX.LT.RIGHT.AND.POSY.LT.TOP.AND. + POSY.GT.BOTTOM) THEN CALL CHRSTR (AXSLAB(J)) ENDIF IF (J.EQ.3) GO TO 62 I1=0 J=J+1 IF (J.EQ.2) THEN I2=INDLAB(2) I3=0 ENDIF IF (J.EQ.3) THEN I3=INDLAB(3) I2=0 ENDIF GO TO 61 62 CONTINUE CALL CHRSIZ (2.5,2.5) CALL DASHLN (0.0,2.0,1.0,0.2) RETURN END SUBROUTINE LCFPREC (ARRAY,IHMX21,IKMX21,ILMAX,NSSQ4) C ============================================== C C---- front end to TJO "precession creator" program C DIMENSION ARRAY(2,IHMX21,IKMX21,0:ILMAX) PARAMETER(MCOLS=500) DIMENSION RSYM(4,4,192) DIMENSION IN(3),ADATA(MCOLS), + DC1(3),DC2(3),DC3(3),AJUNK2(4,4),AJUNK(4,4),IHH(3) C C---- LCFPREC file variables C CHARACTER*80 TITLE CHARACTER*10 PLSTRG CHARACTER*1 CHNNM(24) LOGICAL SQUARE,CALFLG,HPFLG,EOF,LOGMSS(MCOLS) REAL SMIN,SMAX,RLIM,PROJ,PRJLIM INTEGER NSYMB,NTERM,NSPGRP,NMULT,IHPROJ(3) INTEGER ILMAX,ISSQ4 INTEGER HKLMAX,HMIN,HMAX,KMIN,KMAX,LMIN,LMAX,HDIF,LDIF,KDIF INTEGER HHMAX,KKMAX,LLMAX INTEGER H,K,L,DIRECTION,SECTION REAL TOP,BOTTOM,LEFT,RIGHT REAL FMIN,FMAX,AHO,AKO,ALO,FRC,SSQSTP,SSQ4 REAL SMFOBS(1000),SMROBS(1000) INTEGER NRFALL(1000),NRFOBS(1000),LUNIN,LUNOUT C C---- device size 0. VT640 1. HPA4 2. HPA3/CALCOMP 3. HP ACETATE 4. TEKT C DIMENSION DEVCEN(0:4,2),DEVSCL(0:4) CHARACTER HKLZON(20)*3 EXTERNAL MATMULB C C---- Common Blocks C COMMON /SYMMM/ NSYM,NMULT,RSYM COMMON /PINOUT/ LUNIN,LUNOUT COMMON /HKLPLTC/ TITLE,HKLZON COMMON /HKLPLT1/ IFP,SQUARE,CALFLG,HPFLG COMMON /HKLPLT2/ NSYMB,SMIN,SMAX,RLIM,RLIM2,NTERM,NSPGRP, + HKLPDB, FMXX,SIG,IHPROJ,PROJ,PRJLIM,LIST,ANOM,FSCAL,IGEN COMMON /ORTHOG/ RO(4,4),RF(4,4) COMMON /HKLPLT3/RR(3,3,6) C C---- Save Statement C SAVE C C---- Data Statement C C DATA DEVCEN/82.,69.,100.,120. ,120.,70., 130.,120., 120., 90./ DATA DEVCEN/82.,69.,82.,69. ,83.,69., 82.,69., 82., 69./ DATA CHNNM/'A','B','C','D','E','F','G','H','I','J', + 'K','L','M','N','O','P','Q','R','S','T', + 'U','V','W','X'/ C DIRECTION = 0 HMAX = (IHMX21 - 1) /2 KMAX = (IKMX21 - 1) /2 LMAX=ILMAX HMIN=-HMAX KMIN=-KMAX LMIN=0 HKLMAX=HMAX IF (KMAX.GT.HKLMAX)HKLMAX=KMAX IF (LMAX.GT.HKLMAX)HKLMAX=LMAX HKLMIN=-HKLMAX EDGE = 10 SSQSTP=(SMAX -SMIN)/NSSQ4 FMIN=1E20 FMAX=-1E20 HPFLG=.FALSE. CALFLG=.FALSE. C DO I = 1,1000 SMFOBS(I)=0 SMROBS(I)=0 NRFOBS(I)=0 END DO C ISET = -99999999 IF (IGEN .EQ.1) THEN FMAX = 10 FMIN = 0 ISET = 1 NREFLN = (ILMAX+1) * IKMX21 * IHMX21 SUMII = NREFLN SUMSDII = NREFLN WRITE (LUNOUT,'(/,I10)') NREFLN DO 1811 L=0,ILMAX DO 181 K=1,IKMX21 DO 18 IH=1,IHMX21 IHH(1) = IH - HMAX - 1 IHH(2) = K - KMAX - 1 IHH(3) = L C C---- EPS = multiplicity of this reflection. C ISYSAB = 1 - systematic absence C CALL EPSLON(IHH,EPS,ISYSAB) IF (ISYSAB .NE. 1) ARRAY(1,IH,K,L)=ISET IF (ISYSAB .EQ. 1) ARRAY(1,IH,K,L)=-99999999 18 CONTINUE 181 CONTINUE 1811 CONTINUE ENDIF DO 1911 K=0,ILMAX DO 191 J=1,IKMX21 DO 19 I=1,IHMX21 ARRAY(1,I,J,K)=ISET ARRAY(2,I,J,K)=0 19 CONTINUE 191 CONTINUE 1911 CONTINUE C C---- need to read in the data with subroutine below C ok so we now know about the data C display to user and ask what he wants C NAXIS use routine slice C LAYER - use routine slice C 666 CONTINUE IF (IGEN .EQ. 0) THEN HDIF=1 KDIF=1 LDIF=1 IATN = 0 FMIN=999999999.0 FMAX=-9999999999.0 C C---- this reads one piece of data from the MTZ file C and returns it in IDATA - (100 integer array) C 60 CONTINUE CALL LRREFF(1,S,ADATA,EOF) IF(EOF) GO TO 500 CALL LRREFM(1,LOGMSS) NREFL=NREFL+1 IN(1)=NINT(ADATA(1)) IN(2)=NINT(ADATA(2)) IN(3)=NINT(ADATA(3)) IF (S.LT.SMIN .OR. S.GT.SMAX) GO TO 60 C No occurrence of FP IF (LOGMSS(4)) GO TO 60 C SIG is a check whether SIGFP has been set IF (SIG.NE.0.0) THEN IF ( LOGMSS(5)) GO TO 60 IF (ABS(ADATA(5)).LE.0.000000001) GO TO 60 IF( ADATA(4).LT.SIG*ADATA(5)) GO TO 60 SDF=ADATA(5) END IF C IF (.NOT.LOGMSS(6)) THEN IF (ABS(ADATA(6)).GE.0.000000001) 1 ADATA(4) = ADATA(4) + ANOM*ADATA(6) ENDIF NREF2=NREF2+1 FF =ADATA(4) IF (FF.GT.FMAX) THEN FMAX=FF HHMAX=IN(1) KKMAX=IN(2) LLMAX=IN(3) END IF IF (FF.LT.FMIN) FMIN=FF SUMII = SUMII + FF*FF SUMSDII = SUMSDII + 2.0*FF*SDF C C---- so the indices have been found C now stuff the data 's' into the relevant places in the C the array ie I+HMAX+1,K+KMAX+1,L C generate symm equivs C DO 400 J=1,NSYM IH=IN(1)*RSYM(1,1,J) +IN(2)*RSYM(2,1,J) +IN(3)*RSYM(3,1,J) IK=IN(1)*RSYM(1,2,J) +IN(2)*RSYM(2,2,J) +IN(3)*RSYM(3,2,J) IL=IN(1)*RSYM(1,3,J) +IN(2)*RSYM(2,3,J) +IN(3)*RSYM(3,3,J) C C---- always keep L ge 0 C IF (IL.LT.0) THEN IH=-IH IK=-IK IL=-IL ENDIF L=IL IF (IH+HMAX+1.GT.IHMX21 .OR.IH+HMAX+1.LT.1 .OR. + IK+KMAX+1.GT.IKMX21 .OR.IK+KMAX+1.LT.1 .OR. + IL.GT.ILMAX .OR.IL.LT.0) THEN GO TO 411 ENDIF ARRAY(1,IH+HMAX+1,IK+KMAX+1,IL)=FF ARRAY(2,IH+HMAX+1,IK+KMAX+1,IL)=SDF IF ( IL.EQ.0) THEN IH=-IH IK=-IK ARRAY(1,IH+HMAX+1,IK+KMAX+1,IL)=FF ARRAY(2,IH+HMAX+1,IK+KMAX+1,IL)=SDF END IF C HKLPDB.. IF(HKLPDB.GT.0.0)THEN AHO=IH*RF(1,1)+IK*RF(2,1)+IL*RF(3,1) AKO=IH*RF(1,2)+IK*RF(2,2)+IL*RF(3,2) ALO=IH*RF(1,3)+IK*RF(2,3)+IL*RF(3,3) AHO = 100.0*AHO AKO = 100.0*AKO ALO = 100.0*ALO IATN = IATN + 1 OCC = 1.0 /SQRT(S) B = ADATA(4) * 100.0/FMXX WRITE(20,1002) iatn,IH,IK,CHNNM(J),IL,AHO,AKO,ALO,OCC,b 1002 FORMAT('ATOM ',i5,1x,I4,i4,1x,A1,I4,4X,3F8.3,2F6.2,I4) END IF 400 CONTINUE GO TO 60 411 write(LUNOUT,6789)in,ih,ik,il,hmax,kmax,lmax,ff 6789 FORMAT (' H K L outside limits - ?? cell?? ',/, + ' Input H K L- Symmetry H K L - HMAX KMAX LMAX - FF', + /,4X,3I4,4X,3I4,4X,3I4,F20.3) GO TO 60 500 CONTINUE C C---- end of IGEN = 0 C ENDIF C C---- go here if no more data found C WRITE (LUNOUT,'(/,A,2I10)') ' Total number of reflections= ', + NREFL,NREF2 IF (IGEN.EQ.0) THEN RATIO = SUMSDII/SUMII FFMAX=FMAX FFMIN=FMIN WRITE (LUNOUT,'(/,A,2F20.4)') + ' Maximum and minimum values of F = ',FFMAX,FFMIN WRITE (LUNOUT,'(/,A,3i4)') + ' Indices of maximum F = ',HHMAX,KKMAX,LLMAX WRITE (LUNOUT,'(/,A,F10.4)') + ' Ratio of (standard deviation of intensity)/intensity',RATIO WRITE (LUNOUT,'(/,A)') + ' Number of reflections v total possible in recip hemisphere' DO 50211 IL=LMIN,LMAX L=IL DO 5021 IK=KMIN,KMAX K=IK+KMAX +1 DO 502 IH=HMIN,HMAX C Only count +h +-k 0... IF(IL .EQ.0 .AND. IH.LT.0 )GO TO 502 C---- EPS = multiplicity of this reflection. C ISYSAB = 1 - systematic absence C IHH(1) = IH IHH(2) = IK IHH(3) = IL CALL EPSLON(IHH,EPS,ISYSAB) IF (ISYSAB .EQ. 1) GO TO 502 AHO=IH*RF(1,1)+IK*RF(2,1)+IL*RF(3,1) AKO=IH*RF(1,2)+IK*RF(2,2)+IL*RF(3,2) ALO=IH*RF(1,3)+IK*RF(2,3)+IL*RF(3,3) c CALL ORTMAT(RF,IH,IK,IL,AHO,AKO,ALO) SSQ4 = AHO*AHO +AKO*AKO +ALO*ALO - SMIN IF (SSQ4.GT.SMAX) GO TO 502 H=IH+HMAX +1 ISSQ4=INT(SSQ4/SSQSTP) + 1 IF (ISSQ4.GT.NSSQ4)ISSQ4=NSSQ4 IF (ARRAY(1,H,K,L).GT.-1.) THEN SMFOBS(ISSQ4)=SMFOBS(ISSQ4) + ARRAY(1,H,K,L) SUMII = SUMII + FF*FF SUMSDII = SUMSDII + 2.0*FF*SDF SMROBS(ISSQ4)=SMROBS(ISSQ4) + + ARRAY(1,H,K,L)*ARRAY(1,H,K,L)/(2.0*ARRAY(2,H,K,L)*ARRAY(1,H,K,L) + + ARRAY(2,H,K,L)*ARRAY(2,H,K,L)) NRFOBS(ISSQ4)=NRFOBS(ISSQ4) + 1 END IF NRFALL(ISSQ4)=NRFALL(ISSQ4) + 1 502 CONTINUE 5021 CONTINUE 50211 CONTINUE WRITE (LUNOUT,5532) 5532 FORMAT(/1X, $' $TABLE: Data analysis v resln :',/, $' $GRAPHS: Mean Fsq :N:1,3:',/, $' : Mean Fsq/SD :N:1,4:',/, $' : Number of observations, Number possible:N:1,5,6:',/, $' : Completeness v range, Cumulative completeness:N:1,7,8:',/, $' $$'/ +' 1/resol^2 Dmin Mean_FSQ Mean_FSQ/SD Nobs Npos Fract', +' Totfrc $$'/' $$') D0=1000 NRFO=0 NRFA=0 SMFO=0 SMRO=0 SSQ40=SMIN DO 503 IH=1,NSSQ4 SSQ41=SSQ40+SSQSTP IF (SSQ40.GT.0)D0=1./SQRT(SSQ40) IF (SSQ41.GT.0)D1=1./SQRT(SSQ41) FRC=0.0 C Already dealt with multipilicity C NRFALL(IH)=NRFALL(IH)/NMULT SMFO = SMFO + SMFOBS(IH) SMRO = SMRO + SMROBS(IH) NRFA = NRFA + NRFALL(IH) NRFO = NRFO + NRFOBS(IH) IF (NRFALL(IH).GT.0) THEN IF (NRFOBS(IH).GT.0)SMFOBS(IH)=1000.*SMFOBS(IH)/NRFOBS(IH) IF (NRFOBS(IH).GT.0)SMROBS(IH)=SMROBS(IH)/NRFOBS(IH) FRC=FLOAT(NRFOBS(IH))/NRFALL(IH) FRCA=FLOAT(NRFO)/NRFA c WRITE (21,5533)D0,D1,SMFOBS(IH),SMROBS(IH),NRFOBS(IH),NRFALL(IH), c + FRC,FRCA SSQ = 0.5*(SSQ40 +SSQ41) IF (NRFOBS(IH).GT.0) + WRITE (LUNOUT,5533)SSQ,D1,SMFOBS(IH),SMROBS(IH),NRFOBS(IH), + NRFALL(IH),FRC,FRCA 5533 FORMAT (F8.3,F6.2,F12.2,F11.2,2I8,2X,F8.3,4X,F8.3) ENDIF SSQ40=SSQ41 503 CONTINUE WRITE(LUNOUT,*)' $$ ' SMFO = SMFO/NRFO SMRO = SMRO/NRFO WRITE(LUNOUT,5534)SMFO,SMRO,NRFO,NRFA,FRCA 5534 FORMAT (' Totals:',F18.1,F11.2,2I8,14X,F8.3) ENDIF C C---- now I set up my variables C CALL PLOTSETTING (NTERM,TOP,BOTTOM,LEFT,RIGHT,DEVSCL,DEVCEN) ISEC = 1 88 CONTINUE IF (IHPROJ(1).NE.0 .OR. IHPROJ(2).NE.0 .OR. IHPROJ(3).NE.0) + DIRECTION = -1 C C---- use DIRECTION = -1 to mark a projection plot... C WRITE (LUNOUT,'(/,A,I4)') ' Just before slice',DIRECTION IF (DIRECTION.GE.0 .AND. HKLZON(ISEC) .NE.' ') 1 CALL SLICE (SECTION,DIRECTION,ISEC,CALFLG,HPFLG) ISEC = ISEC + 1 IF (DIRECTION.EQ.0) GOTO 666 IF (DIRECTION.EQ.99) GOTO 99 56 CONTINUE C IF (NTERM.NE.4) WRITE (LUNOUT,5566) TITLE C5566 FORMAT (///,A80) CALL MOVETOB (LEFT,TOP-EDGE/2.0) CALL CHRSTR (TITLE) CALL MOVETOB (LEFT,TOP-EDGE) IF (DIRECTION.GT.0) CALL CHRSTR (' Section ') IF (DIRECTION.LT.0) THEN C C---- work out RO for projection C CALL CHRSTR (' Projection ') C C---- get direction cosines of projection axis C AHO=IHPROJ(1)*RF(1,1)+IHPROJ(2)*RF(2,1)+IHPROJ(3)*RF(3,1) AKO=IHPROJ(1)*RF(1,2)+IHPROJ(2)*RF(2,2)+IHPROJ(3)*RF(3,2) ALO=IHPROJ(1)*RF(1,3)+IHPROJ(2)*RF(2,3)+IHPROJ(3)*RF(3,3) c CALL ORTMAT(RF,IHPROJ(1),IHPROJ(2),IHPROJ(3),AHO,AKO,ALO) DDD= SQRT(AHO*AHO + AKO*AKO + ALO*ALO) DC1(1) = AHO/DDD DC1(2) = AKO/DDD DC1(3) = ALO/DDD C C---- find any 2 other perpendicular axes... C PROD = ABS(DC1(1)*DC1(2)*DC1(3)) IF (PROD.GT. 0.0001) THEN DC2(1) = -DC1(2) DC2(2) = DC1(1) DDD= SQRT(DC2(1)*DC2(1) + DC2(2)*DC2(2)) DC2(1) = DC2(1)/DDD DC2(2) = DC2(2)/DDD DC2(3) = 0.0 ELSE DC2(1) = 0.0 DC2(2) = 0.0 DC2(3) = 0.0 IF (DC1(1).LT.0.0001) THEN DC2(1) = 1.0 GO TO 57 ENDIF IF (DC1(2).LT.0.0001) THEN DC2(2) = 1.0 GO TO 57 ENDIF IF (DC1(3).LT.0.0001) THEN DC2(3) = 1.0 ENDIF 57 ENDIF CALL VPROD(DC1,DC2,DC3) IFLG = -1 AJUNK(4,4)=1.0 DO 58 I=1,3 AJUNK(4,I)=0.0 AJUNK(I,4)=0.0 AJUNK(I,1)=DC1(I) AJUNK(I,2)=DC2(I) AJUNK(I,3)=DC3(I) 58 CONTINUE CALL MATMULB(AJUNK2,RF,AJUNK,4,4,4,IFLG) DO 591 I=1,4 DO 59 J=1,4 RF(J,I) = AJUNK2(J,I) 59 CONTINUE 591 CONTINUE GO TO 50 ENDIF IF (DIRECTION.EQ.1) THEN NCODE=1 WRITE(PLSTRG,'(A,I5)') ' h =',SECTION CALL CHRSTR (PLSTRG) C IF (NTERM.NE.4) WRITE (LUNOUT,5567) SECTION C5567 FORMAT (///,' Section H= ',I5) ELSE IF (DIRECTION.EQ.2) THEN NCODE=2 WRITE(PLSTRG,'(A,I5)') ' k =',SECTION CALL CHRSTR (PLSTRG) C IF (NTERM.NE.4) WRITE (LUNOUT,5568) SECTION C5568 FORMAT (///,' Section K= ',I5) ELSE IF (DIRECTION.EQ.3) THEN NCODE=3 WRITE(PLSTRG,'(A,I5)') ' l =',SECTION CALL CHRSTR (PLSTRG) C IF (NTERM.NE.4) WRITE (LUNOUT,5569) SECTION C5569 FORMAT (///,' Section L= ',I5) ENDIF LAYER=SECTION C C---- reset RO and RF appropriate for the reciprocal lattice C section requested C RO(4,4)=1 DO 444 I=1,3 RO(4,I)=0 RO(I,4)=0 DO 44 J=1,3 RO(J,I)=RR(J,I,NCODE) 44 CONTINUE 444 CONTINUE CALL RBRINV(RO,RF) 50 CONTINUE WRITE (LUNOUT,'(A,A,4(/,5X,4F10.5))') + ' *** Orthogonalisation matrix for H K L ', + ' - Sections up HO plotted*** ', + ((RF(J,I),I=1,4),J=1,4) C C---- hello C CALL GRAPHDRAW (ARRAY,IHMX21,IKMX21,ILMAX, + NREFL,HMAX,KMAX,LMAX,HMIN,KMIN,LMIN, + TOP,BOTTOM,LEFT,RIGHT,FMIN,FMAX,SECTION,DIRECTION, + RF,HKLMAX,DEVCEN,DEVSCL) IF (DIRECTION.GE.0) THEN CALL ENDPIC(1) GOTO 88 ENDIF 99 CONTINUE CALL ENDPLT(1) 990 CONTINUE END C C SUBROUTINE MATMULB(BC,B,C,N1,N2,N3,IFLG) C ======================================= C DIMENSION BC(N1,N3),B(N1,N2),C(N2,N3) INTEGER LUNIN,LUNOUT C COMMON /PINOUT/ LUNIN,LUNOUT C SAVE C C---- matrices passed to subroutine C IF (IFLG.LT.0) GO TO 1064 C C---- matrices read in subroutine C IF (IFLG.EQ.0) + WRITE (LUNOUT,*) + ' Input each row for Matrix 1 on a seperate line' IF (IFLG.EQ.1) + WRITE (LUNOUT,*) ' Matrix 1 equals previous product ' DO 6 I=1,N1 IF (IFLG.EQ.0)READ (LUNIN,*)( B(I,J),J=1,N2) WRITE (LUNOUT,100)( B(I,J),J=1,N2) 100 FORMAT (10(10F9.3,/)) 6 CONTINUE WRITE (LUNOUT,*) + ' Input each row for next Matrix on a seperate line' DO 4 I=1,N2 READ (LUNIN,*)( C(I,J),J=1,N3) WRITE (LUNOUT,100)( C(I,J),J=1,N3) 4 CONTINUE 1064 CONTINUE DO 1 I=1,N1 DO 2 J=1,N3 BC(I,J)=0. DO 3 K=1,N2 BC(I,J)=BC(I,J)+B(I,K)*C(K,J) 3 CONTINUE 2 CONTINUE 1 CONTINUE WRITE (LUNOUT,*)' *** Product [M1]*[M2]*[...] ***' DO 5 I=1,N1 WRITE (LUNOUT,100)( BC(I,J),J=1,N3) 5 CONTINUE RETURN END C C SUBROUTINE MOVETOB (X,Y) C ======================= C REAL X,Y C EXTERNAL GSANCD C CALL GSANCD (X,Y) END C C SUBROUTINE PLOTSETTING(NTERM,TOP,BOTTOM,LEFT,RIGHT,DEVSCL,DEVCEN) C ================================================================= C C---- subrountine to set up plotter device C REAL DEVSCL(0:4),DEVCEN(0:4,2),LEFT LOGICAL SQUARE,CALFLG,HPFLG INTEGER LUNIN,LUNOUT C COMMON /PINOUT/ LUNIN,LUNOUT COMMON /HKLPLT1/ IFP,SQUARE,CALFLG,HPFLG C SAVE C CALL BEGPLT ('PLOT') C C---- decide on output device C EDGE = 10 IF (NTERM.EQ.0) THEN TOP =2.*DEVCEN(0,1)-EDGE RIGHT=2.*DEVCEN(0,2)-EDGE CALL DEVICE (NTERM) ELSE IF (NTERM.EQ.3) THEN TOP =2.0*DEVCEN(3,2)-EDGE RIGHT=2.0*DEVCEN(3,1)-EDGE CALL DEVICE (NTERM) ELSE IF (NTERM.EQ.2) THEN CALFLG=.TRUE. TOP =2.*DEVCEN(2,1)-EDGE RIGHT=2.*DEVCEN(2,2)-EDGE CALL DEVICE (NTERM) ELSE IF (NTERM.EQ.1) THEN HPFLG=.TRUE. TOP =2.*DEVCEN(2,1)-EDGE RIGHT=2.*DEVCEN(2,2)-EDGE CALL DEVICE (NTERM) ELSE IF (NTERM.EQ.11) THEN HPFLG=.TRUE. NTERM=1 TOP =2.*DEVCEN(1,1)-EDGE RIGHT=2.*DEVCEN(1,2)-EDGE CALL DEVICE (NTERM) ELSE IF (NTERM.EQ.12) THEN HPFLG=.TRUE. NTERM=1 TOP =2.*DEVCEN(3,1)-EDGE RIGHT=2.*DEVCEN(3,2)-EDGE CALL DEVICE (NTERM) ENDIF BOTTOM=EDGE LEFT=EDGE SCAL1=0.5*(TOP-BOTTOM) SCAL2=0.5*(RIGHT-LEFT) IF (SCAL2.LT.SCAL1)SCAL1=SCAL2 DEVSCL(NTERM)=SCAL1 RETURN END C C SUBROUTINE PLOTUNSET(NTERM,CALFLG,HPFLG) C ======================================== C C---- subroutine to finish ploting C LOGICAL CALFLG,HPFLG INTEGER NTERM,LUNIN,LUNOUT C COMMON /PINOUT/ LUNIN,LUNOUT C SAVE C CALL DEVICE (4) NTERM=1 CALFLG=.FALSE. HPFLG=.FALSE. WRITE (LUNOUT,*)' Plotting to device finished' RETURN END C -------------------------------------------------------------------- C | All graphics calls here. This uses plot9x - PJ Daly 21st May 1990 | C -------------------------------------------------------------------- C C INTEGER FUNCTION SETDEV (STRING) C ================================ C CHARACTER*(*) STRING INTEGER NTERM,NDEVIC,NDIREC C COMMON /DEVATT/ NTERM,NDEVIC,NDIREC C SAVE C NDEVIC = -1 C IF (STRING(1:2).EQ.'ON' .OR. STRING(1:2).EQ.'on') THEN NTERM = 1 NDEVIC = 3 NDIREC = 1 ELSE NTERM = 0 NDEVIC = 1 NDIREC = 0 ENDIF SETDEV = NDEVIC END C SUBROUTINE SLICE(SECTION,DIRECTION,ISEC,CALFLG,HPFLG) C ================================================ C C---- user choice of procession view and of plotting device C INTEGER SECTION,DIRECTION,LEVEL,TEMPSECTION,TEMPDIR CHARACTER TITLE*80,INPUT*3,HKLZON(20)*3 LOGICAL H,K,L,MINUS,CALFLG,HPFLG INTEGER LUNIN,LUNOUT C COMMON /PINOUT/ LUNIN,LUNOUT COMMON /HKLPLTC/ TITLE,HKLZON C SAVE C TEMPSECTION=SECTION TEMPDIR=DIRECTION C C---- write to screen if you are outputting to TEKT C IF (.NOT.CALFLG .AND. .NOT.HPFLG) THEN C CALL TTYON WRITE (LUNOUT,'(A)') ' ** Type RETURN to clear screen' c READ (LUNIN,'(A)') INPUT C CALL TTYOFF CALL MOVETOB (114.0,10.0) ENDIF 2 CONTINUE C2 CALL TTYON C WRITE (LUNOUT,'(A)') ' Which section do you require' C WRITE (LUNOUT,'(A)') ' ? or finish ' C READ (LUNIN,'(A)')INPUT INPUT = HKLZON(ISEC) WRITE (LUNOUT,'(A,A)') ' HKLZON = ',INPUT C CALL TTYOFF CALL CCPUPC(INPUT) IF (INPUT(1:3).EQ.'FIN') THEN DIRECTION=99 GOTO 99 ENDIF SECTION=0 LEVEL=0 H=.FALSE. K=.FALSE. L=.FALSE. MINUS=.FALSE. DO 1 I=1,3 IF (INPUT(I:I).EQ.' ')GO TO 1 IF (INPUT(I:I).EQ.'H') THEN H=.TRUE. GO TO 1 ELSE IF (INPUT(I:I).EQ.'K') THEN K=.TRUE. GO TO 1 ELSE IF (INPUT(I:I).EQ.'L') THEN L=.TRUE. GO TO 1 ELSE IF (INPUT(I:I).EQ.'-') THEN MINUS=.TRUE. GO TO 1 ENDIF READ(INPUT(I:I),300)LEVEL 300 FORMAT (I1) IF (LEVEL.GE.0.AND.LEVEL.LT.10)SECTION=10*SECTION+LEVEL 1 CONTINUE IF (MINUS)SECTION=0-SECTION IF (.NOT. (H))DIRECTION=1 IF (.NOT. (K))DIRECTION=2 IF (.NOT. (L))DIRECTION=3 IF ((.NOT.(H).AND..NOT.(K)).OR.(.NOT.(K).AND..NOT.(L)).OR. + (.NOT.(H).AND..NOT.(L))) THEN WRITE (LUNOUT,*) ' You have not defined a plane ' GOTO 2 ENDIF 99 RETURN END C C SUBROUTINE VPROD(A,B,C) C ======================== C DIMENSION A(3),B(3),C(3) C C---- C will be the vector product of A x B and will be normalised C to unit length C C(1)=A(2)*B(3) - A(3)*B(2) C(2) = A(3)*B(1) - A(1)*B(3) C(3) = A(1)*B(2) - A(2)*B(1) D=SQRT(C(1)*C(1) + C(2)*C(2) + C(3)*C(3) ) IF (D.EQ.0.0) RETURN C(1)=C(1)/D C(2)=C(2)/D C(3)=C(3)/D RETURN END