Title                  :HP45 Calculator program (enhanced) V1mod7
Keywords               :HREM, Wavelength, Scherzer, Resolution, hex, oct, bin
Computer               :anything with F77
Operating System       :anything with F77
Programming Language   :Fortran 77
Hardware Requirements  :CPU, terminal
Authors                :Tony Pitt, Owen Saxton, Graeme Wood, Mike O'Keefe, ...
Correspondence Address :MAOK@CSA2.LBL.GOV
                       :M.A. O'Keefe, NCEM, LBL 72-213, Berkeley, CA 94720

Abstract: This is a program designed to emulate a scientific calculator
operating in reverse polish mode (originally based on a HP45, but now with
many enhancements, including binary, hex, and octal numbers, and electron
microscope functions like wavelength and resolution).   The HELP command
provides one-line explanations of each command.   Commands are -

 +        -        *        /        SQRT     REP      SIN      COS      TAN
 TYPE     DEL      PI       EXP      CLR      CUBE     SQU      REC      SWOP
 SWAP     CHS      POC      ROLL     LOG      ALOG     ASIN     ACOS     ATAN
 LN       OCT      HEX      X        DISP     INT      COMS     BIN      DEC
 INV      DEG      RAD      D>R      R>D      EXIT     WAVL     KV       SCH
 OPT      RSCH     ROPT     LIMD     LIMA     STO      RCL      HELP
 
---------------------------------------------------------------------------

Title                
  :HP45 Calculator program (enhanced) V1mod7
Keywords               :HREM, Wavelength, Scherzer, Resolution, hex, oct, bin
Computer               :anything with F77
Operating System       :anything with F77
Programming Language   :Fortran 77
Hardware Requirements  :CPU, terminal
Authors                :Tony Pitt, Owen Saxton, Graeme Wood, Mike O'Keefe, ...
Correspondence Address :MAOK@CSA2.LBL.GOV
                       :M.A. O'Keefe, NCEM, LBL 72-213, Berkeley, CA 94720

Description: This is a program designed to emulate a scientific calculator
operating in reverse polish mode (originally based on a HP45, but now with
many enhancements, including binary, hex, and octal numbers, and electron
microscope functions like wavelength and resolution).   The HELP command
provides one-line explanations of each command.   Commands are -

 +        -        *        /        SQRT     REP      SIN      COS      TAN
 TYPE     DEL      PI       EXP      CLR      CUBE     SQU      REC      SWOP
 SWAP     CHS      POC      ROLL     LOG      ALOG     ASIN     ACOS     ATAN
 LN       OCT      HEX      X        DISP     INT      COMS     BIN      DEC
 INV      DEG      RAD      D>R      R>D      EXIT     WAVL     KV       SCH
 OPT      RSCH     ROPT     LIMD     LIMA     STO      RCL      HELP

Exact functions of each command are easily read in the Fortran source, which
has extensive inbuilt documentation.  Extra user-defined commands are easily
added to the source at the designated place.  The program stores input numbers
in a stack, and commands operate on the top one or two of the stack members.
Expressions can be up to 80 keystrokes long and are evaluated in reverse-polish
convention.  E.g. the result of the expression 6 5+ is 11, that of 6 5/ is 1.2,
the result of 5 5+5+3/ is 5.

The calculator is entered in "radian-decimal" mode (the prompt expresses this
fact by taking the form "RAD:DEC>" on entry), but can be changed to other
modes by the commands DEG, BIN, OCT, HEX, DEC and RAD.  In radian mode the
commands SIN, COS, TAN assume that all angles are in radians; DEG puts the
calculator into degree mode, and all trig functions then assume that all angles
are measured in degrees; RAD switches the calculator back from degree mode to
radian mode.  The commands BIN, OCT, HEX and DEC are used to switch the
calculator into binary, octal, hexadecimal and decimal modes respectively.
The mode-change commands can be used within arithmetic expressions; e.g. with
the calculator in decimal mode, the command line 23hex2b+dec will place the
decimal number 23 on the stack, convert to hexidecimal mode, place the number
2b on the stack, add the top two numbers, and convert the result to decimal 66.

The commands WAVL, KV, SCH, OPT, RSCH, ROPT, LIMD, LIMA are useful in computing
quantities used in high resolution transmission electron microscopy.  The WAVL
command gives the electron wavelength in Angstroms assuming that the top number
on the stack is the accelerating voltage of the microscope in kV. The command
KV computes the voltage from the wavelength in Angstrom units.

The commands SCH and OPT compute the Scherzer defocus and "optimum" defocus
respectively, assuming that the top two numbers on the stack are the objective
lens spherical aberration coefficient Cs in mm. and the electron wavelength in
Angstroms.  The Scherzer defocus is defined as the square root of (wavelength
times Cs), and the "optimum defocus" as the square root of (wavelength times
Cs times 1.5).  E.g. the command string  400 WAVL 1.1 OPT  is used to compute
the "optimum" defocus for a 400keV microscope with a Cs of 1.1mm.

The command RSCH (Resolution at SCHerzer defocus) gives the resolution in
Angstroms corresponding to the spatial frequency at which the linear-image
CTF reaches cross-over at Scherzer defocus; the command ROPT (Resolution at
OPTimum defocus) gives the resolution in Angstroms at the spatial frequency
at which the linear-image CTF reaches the upper limit of the optimum defocus
pass-band.  The command string  400 WAVL 1.1 ROPT  will give the "structure-
image" resolution for a 400keV electron microscope with a Cs of 1.1mm.

The commands LIMD and LIMA give the "Information Limit" resolutions due to the
effects of partial temporal and spatial coherence respectively -- i.e. the
limits at which the "envelope functions" due to chromatic aberration (spread
of focus) and incident beam convergence drop to exp(-2).

LIMD is the resolution LIMit (in A) due to Delta, where delta is the spread of
focus in Angstrom units due to the combined effects of chromatic aberration,
high-voltage ripple, energy spread in the electron beam, lens current ripple,
and vertical vibration of the specimen relative to the objective lens polepiece.
LIMD is calculated from the electron wavelength (in A) and the spread of focus
in Angstrom: e.g. the expression  400 WAVL 85 LIMD  would give the spread-of-
focus resolution limit of a 400keV microscope with a delta of 85A as 1.48A.

LIMA is the resolution LIMit (in A) due to Alpha, where alpha is the semi-angle
measured in milliradian of the cone of electrons incident on the specimen.  The
value of LIMA varies with defocus and Cs, and thus requires the presence of
four parameters on the stack -- the electron wavelength (in A), Cs(mm), the
defocus(A), and alpha(milliradian).  The expression  400 WAVL 1.1 -500 .5 LIMA
calculates a convergence resolution limit of 1.27A for a 400keV microscope with
a Cs of 1.1mm focussed at an underfocus of 500A with a convergence semi-angle
of 0.5milliradian.  NOTE: Although the other three parameters may be entered in
arbitrary order, the last entry on the stack must be the value of alpha.

References: O'Keefe & Pitt, Electron Microscopy 1980, vol.1, pp.122-123.

Compilation procedure: standard F77
Linking procedure: standard
Test data: none
General comments: i/o is to Fortran unit 3, so assign this logical unit
                  to your terminal
Title                  :HP45 Calculator program (enhanced) V1mod7
Keywords               :HREM, Wavelength, Scherzer, Resolution, hex, oct, bin
Computer               :anything with F77
Operating System       :anything with F77
Programming Language   :Fortran 77
Hardware Requirements  :CPU, terminal
Authors                :Tony Pitt, Owen Saxton, Graeme Wood, Mike O'Keefe, ...
Correspondence Address :MAOK@CSA2.LBL.GOV
                       :M.A. O'Keefe, NCEM, LBL 72-213, Berkeley, CA 94720


C HP calculator program by TJPitt.
C
C 91.05.06   MAOK - Add functions LIMD and LIMA
C
C 88.01.04   MAOK - Add HELP command to explain commands
C
C 87.10.01   MAOK - Add facilities to STO and RCL to/from 32 registers
C
C 86.10.10   MAOK - Include functions Wavl, kV, Sch, Opt, RSch, ROpt
C
C 85.05.15   GJW - Modified to recognise lower-case commands
C
C 80.10.30   WOS - V1mod1 bug associated with reading an octal number
C            after a decimal number containing an exponent cured
C
C 80.05.09   WOS - Public release version 0
C            Degree/radian handling inserted, and number output
C            routine tidied up.  Minus operator between two numbers
C            (separated from second by space) now handled correctly.
C            Negative exponents now allowable!
C
C 80.02.15   TJP - Started
C
C
C      NCOMM is the number of commands
       PARAMETER (NCOMM=53)
       DOUBLE PRECISION INTEGER KNUM
       DOUBLE PRECISION STACK,RNUM,DR,DR1,DR2,DR3,PI,CONV1,CONV2,R
       DOUBLE PRECISION STORE(32),ALIMIT
       LOGICAL RDNUM,FMINUS,FREGIS,WKDEG
       CHARACTER*26 ALPHAU,ALPHAL
       CHARACTER*3 IFLAG,JFLAG
       CHARACTER*1 INCHAR
       INTEGER IN(80),IDIGIT(16),COMMS(4,NCOMM)
C
       COMMON /HPCALC/ Pi,IBASE,IDIGIT,STACK(64)
C
       DATA ALPHAL/'abcdefghijklmnopqrstuvwxyz'/
       DATA ALPHAU/'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
       DATA IN/'C','O','M','S',76*' '/
       DATA ISTLN/64/,ISPACE/' '/,MINUS/'-'/,PI/3.14159265358979D0/
C CONV1 is multiplier to get from degrees to radians
C CONV2 is the multiplier for the converse
       DATA CONV1/.174532925199433D-1/,CONV2/.572957795130824D2/
C Store 16 Hexadecimal digits
       DATA IDIGIT/'0','1','2','3','4','5','6','7','8','9',
     1  'A','B','C','D','E','F'/
C Here are the commands
       DATA COMMS/'+',3*' ','-',3*' ','*',3*' ','/',3*' ',
     1  'S','Q','R','T','R','E','P',' ','S','I','N',' ','C','O','S',' ',
     2  'T','A','N',' ','T','Y','P','E','D','E','L',' ','P','I',' ',' ',
     3  'E','X','P',' ','C','L','R',' ','C','U','B','E','S','Q','U',' ',
     4  'R','E','C',' ','S','W','O','P','S','W','A','P','C','H','S',' ',
     5  'P','O','C',' ','R','O','L','L','L','O','G',' ','A','L','O','G',
     6  'A','S','I','N','A','C','O','S','A','T','A','N','L','N',' ',' ',
     7  'O','C','T',' ','H','E','X',' ','X',' ',' ',' ','D','I','S','P',
     8  'I','N','T',' ','C','O','M','S','B','I','N',' ','D','E','C',' ',
     9  'I','N','V',' ','D','E','G',' ','R','A','D',' ','D','>','R',' ',
     1  'R','>','D',' ','E','X','I','T','W','A','V','L','K','V',' ',' ',
     2  'S','C','H',' ','O','P','T',' ','R','S','C','H','R','O','P','T',
     3  'L','I','M','D','L','I','M','A',
     4  'S','T','O',' ','R','C','L',' ','H','E','L','P'/
C
C Initialise everything
       DO 1 IS=1,ISTLN
1      STACK(IS)=0.
       DO 14 I=5,80
14     IN(I)=ISPACE
       IBASE=10
       WKDEG=.FALSE.
       IFLAG='Rad'
       JFLAG='Dec'
       ISTPTR=0
       FMINUS=.FALSE.
       FREGIS=.FALSE.
       DO 10 I=1,6
         IDIGIT(I+10)=ICHAR(ALPHAU(I:I))
10     CONTINUE
       DO 13 L=1,NCOMM
         DO 12 K=1,4
           INCHAR=CHAR(COMMS(K,L))
           DO 11 I=1,26
             IF(INCHAR.EQ.ALPHAU(I:I)) THEN
               COMMS(K,L)=ICHAR(ALPHAU(I:I))
               GOTO 12
             END IF
11         CONTINUE
12       CONTINUE
13     CONTINUE
C
C Put out the header and clear unit 0
       WRITE (3,1009)
1009   FORMAT (27X,' ***  HP45 - V1MOD7  ***'/)
C
C Write out the 'COMS' result
       FREGIS=.FALSE.
       GOTO 15
C
C Terminal input here
10000  FREGIS=.FALSE.
       WRITE (3,20) IFLAG,JFLAG
20     FORMAT (' ',A3,':',A3,'> ',$)
       READ (3,21,ERR=10000,END=9999) IN
21     FORMAT (80A1)
15     INPTR=1
       DO 2 I=1,80
       LENGTH=81-I
       IF (IN(LENGTH).NE.ISPACE) GOTO 18
2      CONTINUE
       GOTO 10000
C
C Convert all letters to uppercase
18     CONTINUE
       DO 19 L=1,LENGTH
         INCHAR=CHAR(IN(L))
         DO 17 I=1,26
           IF(INCHAR.EQ.ALPHAL(I:I).OR.INCHAR.EQ.ALPHAU(I:I)) THEN
             IN(L)=ICHAR(ALPHAU(I:I))
             GOTO 19
           END IF
17       CONTINUE
19     CONTINUE
C
C Back to here for the next thing out of the input buffer
C LENGTH points to the last non-space
C INPTR is where to start looking now
C First we must skip spaces
3      IF (IN(INPTR).NE.ISPACE) GOTO 4
       INPTR=INPTR+1
       IF (INPTR.GT.LENGTH) GOTO 10000
       GOTO 3
C
C Next see if we have a '-'
4      FMINUS=IN(INPTR).EQ.MINUS
       IF (.NOT.FMINUS) GOTO 5
       INPTR=INPTR+1
       IF (IN(INPTR).EQ.ISPACE) GOTO 2002
       IF (INPTR.LE.LENGTH) GOTO 5
       GOTO 2002
C
C Try recognising a command
5      IF (FREGIS) GOTO 9
       ICH=IN(INPTR)
       DO 7 IC=1,NCOMM
       IF (ICH.NE.COMMS(1,IC)) GOTO 7
       JPTR=INPTR
       DO 6 I=2,4
       JPTR=JPTR+1
       IF (COMMS(I,IC).EQ.ISPACE) GOTO 8
       IF (JPTR.GT.LENGTH) GOTO 7
       IF (COMMS(I,IC).NE.IN(JPTR)) GOTO 7
6      CONTINUE
       JPTR=JPTR+1
       GOTO 8
7      CONTINUE
C
C Not a command - try a number
9      IF (.NOT.RDNUM(IN,INPTR,RNUM)) GOTO 9002
C Got one - stack it, if there's room
       ISTPTR=ISTPTR+1
       IF (ISTPTR.GT.ISTLN) GOTO 9001
       IF (FMINUS) RNUM=-RNUM
       STACK(ISTPTR)=RNUM
       IF (FREGIS) GOTO 30129
       IF (INPTR.GT.LENGTH) GOTO 10000
       GOTO 3
C
C Match found - IC is pointer to command
8      INPTR=JPTR
       IF (FMINUS) GOTO 2002
30129  GOTO (2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,
     1       2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,
     2       2021,2022,2023,2024,2025,2026,2027,2028,2029,2030,
     3       2031,2032,2033,2034,2035,2036,2037,2038,2039,2040,
     4       2041,2031,2042,2043,2044,2045,2046,2047,2048,2049,
     5       2050,2051,2052),IC
C
C Here are the commands
C
C '+':  add the top two numbers on the stack
2001   IF (ISTPTR.LT.2) GOTO 9003
       ISTPTR=ISTPTR-1
       STACK(ISTPTR)=STACK(ISTPTR)+STACK(ISTPTR+1)
       GOTO 2099
C
C '-':  subtract top two numbers on stack
2002   IF (ISTPTR.LT.2) GOTO 9003
       ISTPTR=ISTPTR-1
       STACK(ISTPTR)=STACK(ISTPTR)-STACK(ISTPTR+1)
       FMINUS=.FALSE.
       GOTO 2099
C
C '*':  multiply top two numbers on stack
2003   IF (ISTPTR.LT.2) GOTO 9003
       ISTPTR=ISTPTR-1
       STACK(ISTPTR)=STACK(ISTPTR)*STACK(ISTPTR+1)
       GOTO 2099
C
C '/':  divide top two numbers on the stack
2004   IF (ISTPTR.LT.2) GOTO 9003
       IF (STACK(ISTPTR)) 20041,9005,20041
20041  ISTPTR=ISTPTR-1
       STACK(ISTPTR)=STACK(ISTPTR)/STACK(ISTPTR+1)
       GOTO 2099
C
C 'SQRT':  square root top number on stack
2005   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20051,20051
20051  STACK(ISTPTR)=DSQRT(STACK(ISTPTR))
       GOTO 2099
C
C 'REP':  repeat top number on stack
2006   IF (ISTPTR.LT.1) GOTO 9003
       ISTPTR=ISTPTR+1
       IF (ISTPTR.GT.ISTLN) GOTO 9001
       STACK(ISTPTR)=STACK(ISTPTR-1)
       GOTO 2099
C
C 'SIN':  take sin of top number on stack
2007   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (WKDEG) RNUM=RNUM*CONV1
       STACK(ISTPTR)=DSIN(RNUM)
       GOTO 2099
C
C 'COS':  take cosine of top number on stack
2008   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (WKDEG) RNUM=RNUM*CONV1
       STACK(ISTPTR)=DCOS(RNUM)
       GOTO 2099
C
C 'TAN':  take tangent of top number on stack
2009   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (WKDEG) RNUM=RNUM*CONV1
       STACK(ISTPTR)=DTAN(RNUM)
       GOTO 2099
C
C 'TYPE':  type the top number on the stack
2010   IF (INPTR.LT.LENGTH) CALL OUTNUM(STACK(ISTPTR))
       GOTO 2099
C
C 'DEL': delete top number from stack
2011   IF (ISTPTR.LT.1) GOTO 9003
       ISTPTR=ISTPTR-1
       GOTO 2099
C
C 'PI': put PI on the stack
2012   ISTPTR=ISTPTR+1
       IF (ISTPTR.GT.ISTLN) GOTO 9001
       STACK(ISTPTR)=PI
       GOTO 2099
C
C 'EXP': raise E to the power of the last number on stack
2013   IF (ISTPTR.LT.1) GOTO 9003
       R=STACK(ISTPTR)
       IF (R.GT.170) GOTO 9006
       IF (R.LT.-170) GOTO 9007
       STACK(ISTPTR)=DEXP(STACK(ISTPTR))
       GOTO 2099
C
C 'CLR': clear the stack
2014   ISTPTR=0
       GOTO 2099
C
C 'CUBE': cube the top number on the stack
2015   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       STACK(ISTPTR)=RNUM*RNUM*RNUM
       GOTO 2099
C
C 'SQU': square the top number on the stack
2016   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       STACK(ISTPTR)=RNUM*RNUM
       GOTO 2099
C
C 'REC': take the reciprocal of the top number on the stack
2017   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 20171,9005,20171
20171  STACK(ISTPTR)=1.D0/RNUM
       GOTO 2099
C
C 'SWOP': see below! - that's because i can't spell swap!
2018   CONTINUE
C
C 'SWAP': exchange the top two numbers on the stack
2019   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       STACK(ISTPTR)=STACK(ISTPTR-1)
       STACK(ISTPTR-1)=RNUM
       GOTO 2099
C
C 'CHS': change sign of top number on stack
2020   IF (ISTPTR.LT.1) GOTO 9003
       STACK(ISTPTR)=-STACK(ISTPTR)
       GOTO 2099
C
C 'POC': raise next to top number to the power of the top one
2021   IF (ISTPTR.LT.2) GOTO 9003
       DR1=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       DR=STACK(ISTPTR)
       IF(DR.EQ.0) GOTO 20212
       IF (DR.LT.0.D0) GOTO 20211
       STACK(ISTPTR)=DEXP(DLOG(DR)*DR1)
       GOTO 2099
20211  IR1=INT(DR1)
       IF((DR1-IR1).NE.0) GOTO 9010
       STACK(ISTPTR)=DEXP(DLOG(-DR)*DR1)
       IF((IR1-2*(IR1/2)).NE.0) STACK(ISTPTR)=-STACK(ISTPTR)
       GOTO 2099
20212  STACK(ISTPTR)=0.D0
       GOTO 2099
C
C 'ROLL': roll the stack, pushing top number to the bottom
2022   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       J=ISTPTR
20221  STACK(J)=STACK(J-1)
       J=J-1
       IF (J.GT.1) GOTO 20221
       STACK(1)=RNUM
       GOTO 2099
C
C 'LOG': take the logarithm of top number on the stack
2023   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,9009,20231
20231  STACK(ISTPTR)=DLOG10(RNUM)
       GOTO 2099
C
C 'ALOG': take anti-logarithm of top number on stack
2024   IF (ISTPTR.LT.1) GOTO 9003
       R=STACK(ISTPTR)
       IF (R.GT.70D0) GOTO 9006
       IF (R.LT.-70D0) GOTO 9007
       STACK(ISTPTR)=10.**R
       GOTO 2099
C
C 'ASIN': arcsine of top number on stack
2025   IF (ISTPTR.LT.1) GOTO 9003
       R=STACK(ISTPTR)
       IF (DABS(R).GT.1.) GOTO 9008
       RNUM=DASIN(R)
       IF (WKDEG) RNUM=RNUM*CONV2
       STACK(ISTPTR)=RNUM
       GOTO 2099
C
C 'ACOS': arccosine of top number on stack
2026   IF (ISTPTR.LT.1) GOTO 9003
       R=STACK(ISTPTR)
       IF (ABS(R).GT.1D0) GOTO 9008
       RNUM=DACOS(R)
       IF (WKDEG) RNUM=RNUM*CONV2
       STACK(ISTPTR)=RNUM
       GOTO 2099
C
C 'ATAN': arctangent of top number on stack
2027   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=DATAN(STACK(ISTPTR))
       IF (WKDEG) RNUM=RNUM*CONV2
       STACK(ISTPTR)=RNUM
       GOTO 2099
C
C 'LN': natural logarithm of top number on stack
2028   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,9009,20281
20281  STACK(ISTPTR)=DLOG(RNUM)
       GOTO 2099
C
C 'OCT': i/o to octal
2029   IBASE=8
       JFLAG='Oct'
       GOTO 2099
C
C 'HEX': i/o to hexadecimal
2030   IBASE=16
       JFLAG='Hex'
       GOTO 2099
C
C 'EXIT': we're done - quit
2031   WRITE(3,20311)
20311  FORMAT(30X,'Exit')
       GOTO 9999
C
C 'DISP': type the whole stack
2032   IF (ISTPTR.LE.1) GOTO 2099
       ITMP01=ISTPTR-1
       DO 20321 J=1,ITMP01
20321  CALL OUTNUM (STACK(J))
       GOTO 2099
C
C 'INT': take integer part of top number on stack
2033   IF (ISTPTR.LT.1) GOTO 9003
       KNUM=IDINT(STACK(ISTPTR))
       STACK(ISTPTR)=KNUM
       GOTO 2099
C
C 'COMS': type out the list of commands
2034   WRITE (3,20341) ((COMMS(I,J),I=1,4),J=1,NCOMM)
20341  FORMAT (7(3X,4A1,4X))
       WRITE (3,*) ' '
       GOTO 3
C
C 'BIN': i/o in binary directly
2035   IBASE=2
       JFLAG='Bin'
       GOTO 2099
C
C 'DEC': restore decimal i/o mode
2036   IBASE=10
       JFLAG='Dec'
       GOTO 2099
C
C 'INV':  same as 'REC' - reciprocal
2037   GOTO 2017
C
C 'DEG': convert to degrees mode
2038   WKDEG=.TRUE.
       IFLAG='Deg'
       GOTO 3
C
C 'RAD': convert to radian mode
2039   WKDEG=.FALSE.
       IFLAG='Rad'
       GOTO 3
C
C 'D>R': convert degrees to radians
2040   IF (ISTPTR.LT.1) GOTO 9003
       STACK(ISTPTR)=STACK(ISTPTR)*CONV1
       GOTO 2099
C
C 'R>D': convert radians to degrees
2041   IF (ISTPTR.LT.1) GOTO 9003
       STACK(ISTPTR)=STACK(ISTPTR)*CONV2
       GOTO 2099
C
C 'WAVL': compute electron wavelength in A from kV
2042   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20421,20421
20421  DR=STACK(ISTPTR)
       STACK(ISTPTR)=0.387818D0/DSQRT(DR*(1.D0+.978459D-3*DR))
       GOTO 2099
C
C 'kV': compute microscope voltage in kV from wavelength in A
2043   IF (ISTPTR.LT.1) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,9009,20431
20431  DR=STACK(ISTPTR)
       STACK(ISTPTR) = 511.007615D0
     &                 *(DSQRT(1.0D0+(0.02426214948D0/DR)**2)-1.0D0)
       GOTO 2099
C
C 'SCH': Scherzer defocus from wavl(A) and Cs(mm) (top 2 stack values)
2044   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20441,20441
20441  DR=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20442,20442
20442  DR1=STACK(ISTPTR)
       STACK(ISTPTR)=-DSQRT(1.D7*DR*DR1)
       GOTO 2099
C
C 'OPT': Optimum defocus (= sqrt(3/2) times Scherzer defocus)
2045   STACK(ISTPTR)=1.5D0*STACK(ISTPTR)
       GOTO 2044
C
C 'RSCH': Resolution at Scherzer defocus (where sin(chi) crosses axis)
2046   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20461,20461
20461  DR=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20462,20462
20462  DR1=STACK(ISTPTR)
       IF (DR.GE.DR1) Then
         STACK(ISTPTR) = DSQRT(DSQRT(DR*DR1*DR1*DR1*1.0D7)/2.0D0)
       Else
         STACK(ISTPTR) = DSQRT(DSQRT(DR1*DR*DR*DR*1.0D7)/2.0D0)
       EndIf
       GOTO 2099
C
C 'ROPT': Resolution at Optimum defocus (where sin(chi) is -.707)
2047   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20471,20471
20471  DR=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20472,20472
20472  DR1=STACK(ISTPTR)
       IF (DR.GE.DR1) Then
         STACK(ISTPTR) = DSQRT(DSQRT(DR*DR1*DR1*DR1*1.0D7))/1.4916D0
       Else
         STACK(ISTPTR) = DSQRT(DSQRT(DR1*DR*DR*DR*1.0D7))/1.4916D0
       EndIf
       GOTO 2099
C
C 'LIMD': resolution LIMit imposed by Delta
C         i.e. where the spread-of-focus envelope falls to exp(-2)
C  Ref: O'Keefe & Pitt, Electron Microscopy 1980, vol.1, pp.122-123.
2048   IF (ISTPTR.LT.2) GOTO 9003
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20481,20481
20481  DR=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       RNUM=STACK(ISTPTR)
       IF (RNUM) 9009,20482,20482
20482  DR1=STACK(ISTPTR)
       STACK(ISTPTR) = dsqrt(0.5D0*PI*DR*DR1)
       GOTO 2099

C 'LIMA': resolution LIMit imposed by Alpha
C         i.e. where the convergence envelope falls to exp(-2)
2049   IF (ISTPTR.LT.4) GOTO 9003
       DR=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       DR1=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       DR2=STACK(ISTPTR)
       ISTPTR=ISTPTR-1
       DR3=STACK(ISTPTR)
       STACK(ISTPTR) = ALIMIT(DR,DR1,DR2,DR3)
       RNUM = STACK(ISTPTR)
       IF(RNUM) 9008,9008,2099
C
C 'STO n': Store current value in STORE(n)
2050   IF (FREGIS) THEN
          FREGIS=.FALSE.
          IF (ISTPTR.LT.2) GOTO 9003
          IREGIS=INT(STACK(ISTPTR))+1
          ISTPTR=ISTPTR-1
          IF (IREGIS.LT.1.OR.IREGIS.GT.32) GOTO 9011
          STORE(IREGIS)=STACK(ISTPTR)
       ELSE
          FREGIS=.TRUE.
       END IF
       GOTO 2099
C
C 'RCL n': Recall a value from STORE(n)
2051   IF (FREGIS) THEN
          FREGIS=.FALSE.
          IF (ISTPTR.LT.1) GOTO 9003
          IREGIS=INT(STACK(ISTPTR))+1
          ISTPTR=ISTPTR-1
          IF (IREGIS.LT.1.OR.IREGIS.GT.32) GOTO 9011
          ISTPTR=ISTPTR+1
          STACK(ISTPTR)=STORE(IREGIS)
       ELSE
          FREGIS=.TRUE.
       END IF
       GOTO 2099
C
C  'HELP': Write out a list of all commands
2052   WRITE (3,*) ' '
       WRITE (3,20341) ((COMMS(I,J),I=1,4),J=1,NCOMM)
       WRITE (3,*) ' '
       WRITE (3,20342)
20342  FORMAT (
     $ 6x,'+   :  add the top two numbers on the stack',/
     $,6x,'-   :  subtract top two numbers on the stack',/
     $,6x,'*   :  multiply top two numbers on stack',/
     $,6x,'/   :  divide top two numbers on the stack',/
     $,6x,'SQRT:  square root top number on stack',/
     $,6x,'REP :  repeat top number on stack',/
     $,6x,'SIN :  take sin of top number on stack',/
     $,6x,'COS :  take cosine of top number on stack',/
     $,6x,'TAN :  take tangent of top number on stack',/
     $,6x,'TYPE:  type the top number on the stack')
       WRITE (3,20350)
20350  FORMAT (/,' (more....)')
       READ (3,21) INCHAR
       WRITE (3,20343)
20343  FORMAT (
     $ 6x,'DEL :  delete top number from stack',/
     $,6x,'PI  :  put PI on the stack',/
     $,6x,'EXP :  raise E to the power of the last number on stack',/
     $,6x,'CLR :  clear the stack',/
     $,6x,'CUBE:  cube the top number on the stack',/
     $,6x,'SQU :  square the top number on the stack',/
     $,6x,'REC :  take the reciprocal of the top number on the stack',/
     $,6x,'SWOP:  see below! - that''s because i can''t spell swap!',/
     $,6x,'SWAP:  exchange the top two numbers on the stack',/
     $,6x,'CHS :  change sign of top number on stack',/
     $,6x,'POC :  raise next to top number to the power '
     1,   'of the top one',/
     $,6x,'ROLL:  roll the stack, pushing top number to the bottom',/
     $,6x,'LOG :  take the logarithm of top number on the stack',/
     $ 6x,'ALOG:  take anti-logarithm of top number on stack',/
     $,6x,'ASIN:  arcsine of top number on stack',/
     $,6x,'ACOS:  arccosine of top number on stack',/
     $,6x,'ATAN:  arctangent of top number on stack',/
     $ 6x,'LN  :  natural logarithm of top number on stack',/
     $,6x,'OCT :  i/o to octal',/
     $,6x,'HEX :  i/o to hexadecimal',/
     $,6x,'EXIT:  we''re done - so quit the program')
       WRITE (3,20350)
       READ (3,21) INCHAR
       WRITE (3,20345)
20345  FORMAT (
     $ 6x,'DISP:  type the whole stack',/
     $,6x,'INT :  take integer part of top number on stack',/
     $,6x,'COMS:  list the commands',/
     $,6x,'BIN :  i/o in binary directly',/
     $,6x,'DEC :  restore decimal i/o mode',/
     $,6x,'INV :  same as ''REC'' - reciprocal',/
     $,6x,'DEG :  convert to degrees mode',/
     $,6x,'RAD :  convert to radian mode',/
     $,6x,'D>R :  convert degrees to radians',/
     $,6x,'R>D :  convert radians to degrees',/
     $,6x,'STO n: Store current value in STORE(n)',/
     $,6x,'RCL n: Recall a value from STORE(n)',/
     $,6x,'HELP:  type out this list of commands')
       WRITE (3,20350)
       READ (3,21) INCHAR
       WRITE (3,20346)
20346  FORMAT (
     $ 6x,'WAVL:  compute electron wavelength in A from kV',/
     $,6x,'kV  :  compute microscope voltage in kV '
     1,   'from wavelength in A',/
     $,6x,'SCH :  Scherzer defocus from top two stack values '
     1,   'as wavl(A), Cs(mm)',/
     $,6x,'OPT :  Optimum defocus from top two stack values '
     1,   'as wavl(A), Cs(mm)',/,13x,'(= sqrt(3/2) times '
     2,   'Scherzer defocus)',/
     $,6x,'RSCH:  Resolution at Scherzer from top two stack values '
     1,   'as wavl(A), Cs(mm)',/,13x,'(where sin(chi) crosses axis)',/
     $,6x,'ROPT:  Resolution at Optimum from top two stack values '
     1,   'as wavl(A), Cs(mm)',/,13x,'(where sin(chi) is -.707)',/
     $,6x,'LIMD:  resolution LIMit (in A) due to Delta -- '
     1,   'from top two stack values',/,13x,'as wavl(A) and delta(A) '
     2,   '-- (where the spread-of-focus envelope',/,13x,'falls to a '
     3,   'value of exp(-2)) ',/
     $,6x,'LIMA:  resolution LIMit (in A) due to Alpha -- '
     1,   'from top four stack values',/,13x,'as wavl(A), Cs(mm), '
     2,   'defocus(A), and alpha(milliradian).  Alpha is ',/,13x
     3,   'the half-angle of the convergence disk, and '
     4,   'MUST be entered last.',/,13x,'-- (where the convergence '
     5,   'envelope falls to a value of exp(-2))')
       WRITE (3,*) ' '
       WRITE (3,*) ' '
       WRITE (3,*) ' '
       GOTO 3
C
C ****** ****** ****** ****** ****** ****** ******
C
C Insert more commands here when the others work
C
C ****** ****** ****** ****** ****** ****** ******
C
C Here is the end of a command, where we type out the
C Top number on the stack - usually the result of
C The last command
2099   IF (FMINUS) GOTO 30129
       IF (INPTR.LE.LENGTH) GOTO 3
       IF (ISTPTR.EQ.0) GOTO 9004
       CALL OUTNUM (STACK(ISTPTR))
       GOTO 10000
C
C Here are the error messages - all assume subsequently
C That we are going to take a new line of input
C
C Stack full
9001   WRITE(3,9101)
9101   FORMAT(30X,'Stack full')
       ISTPTR=ISTPTR-1
       GOTO 10000
C
C Syntax
9002   WRITE(3,9102)
9102   FORMAT(30X,'Syntax')
       GOTO 10000
C
C Insufficient values on the stack
9003   WRITE(3,9103)
9103   FORMAT(30X,'Insufficient values on stack')
       GOTO 10000
C
C Stack empty
9004   WRITE(3,9104)
9104   FORMAT(30X,'Stack empty')
        GOTO 3
C
C Divide by zero
9005   WRITE(3,9105)
9105   FORMAT(30X,'Divide by zero')
       GOTO 10000
C
C Overflow
9006   WRITE(3,9106)
9106   FORMAT(30X,'Overflow')
       GOTO 10000
C
C Underflow
9007   WRITE(3,9107)
9107   FORMAT(30X,'Underflow')
       STACK(ISTPTR)=0.
       CALL OUTNUM(0.)
       GOTO 10000
C
C Outside range
9008   WRITE(3,9108)
9108   FORMAT(30X,'Outside range')
       GOTO 10000
C
C Less than zero
9009   WRITE(3,9109)
9109   FORMAT(30X,'Less than or equal to zero')
       GOTO 10000
C
C -ve number to power of non-integer
9010   WRITE(3,9110)
9110   FORMAT(30X,'Negative number to power of non-integer')
       ISTPTR=ISTPTR+1
       GOTO 10000
C
C Store or Recall outside range
9011   WRITE(3,9111)
9111   FORMAT(30X,'Only registers 0 through 31 are valid')
       GOTO 10000
C
C End or err on input - just quit tidily!
9999   CONTINUE
C
       END
C RDNUM - Read a positive number in various ways
C
       LOGICAL FUNCTION RDNUM (IN,INPTR,RNUM)
C
C IN Input buffer
C INPTR  Pointer in the buffer, advanced past number
C RNUM The number read
C RDNUM returns TRUE if a number was found, false otherwise
C
C Additional mods at Cambridge to make it read fractions
C and exponents in any base
C
       DOUBLE PRECISION RNUM,R,FMUL1,RBASE,Pi
       LOGICAL FNUM,DEC,FRAC,EXPON,FMINUS
       INTEGER IN(80),IDIGIT(16)
C
       COMMON /HPCALC/ Pi,IBASE,IDIGIT,STACK(64)
C
       EQUIVALENCE (IE,IDIGIT(15))
C
       DATA IMINUS/'-'/,IPT/'.'/
C
       FRAC=.FALSE.
       FMINUS=.FALSE.
       FMUL1=1.
       FNUM=.FALSE.
      EXPON=.FALSE.
       R=0.
      DEC=IBASE.EQ.10
       RBASE=DFLOAT(IBASE)
C
C Reading in octal or hexadecimal or binary
801    J=IN(INPTR)
805    DO 81 ID=1,IBASE
       IF (J.EQ.IDIGIT(ID)) GOTO 82
81     CONTINUE
       GOTO 90
C Found a legal digit
82     FNUM=.TRUE.
       R=R*IBASE+ID-1.
       IF (FRAC)FMUL1=FMUL1/IBASE
       INPTR=INPTR+1
       IF (INPTR.GT.80) GOTO 90
       GOTO 801
C
C Finished reading our number
90     IF (FRAC) R=R*FMUL1
      GOTO 101
900    IF (FNUM) GOTO 902
C See if we found a number
901    RDNUM=.FALSE.
       RETURN
C
C And we found one
902    RNUM=R
104    RDNUM=.TRUE.
       RETURN
C
C Come to the end (for the moment) of our number
C but first we need to check for various things
101    IF (INPTR.GT.80) GOTO 900
       IF (J.NE.IPT) GOTO 102
       IF (EXPON.OR.FRAC) GOTO 103
       INPTR=INPTR+1
       IF (INPTR.GT.80) GOTO 103
       FRAC=.TRUE.
       GOTO 801
C Not a point - try an exponent
102    IF (.NOT.FNUM) GOTO 901
       IF (J.NE.IE.OR.EXPON) GOTO 103
       RNUM=R
       R=0.
       FRAC=.FALSE.
       FNUM=.FALSE.
       EXPON=.TRUE.
       INPTR=INPTR+1
       IF (INPTR.GT.80) GOTO 106
       J=IN(INPTR)
       FMINUS=J.EQ.IMINUS
       IF (.NOT.FMINUS) GOTO 805
       INPTR=INPTR+1
       IF (INPTR.LE.80) GOTO 801
       INPTR=INPTR-1
106    INPTR=INPTR-1
       GOTO 104
C
C Okay we're done - return
103    IF (.NOT.EXPON) GOTO 902
       IF (.NOT.FNUM) GOTO 106
       I=R
       IF (FMINUS) I=-I
       RNUM=RNUM*(rbase**i)
       GOTO 104
C
       END
C REPLACEMENT OUTNUM FOR TONY'S HP CALCULATOR PROGRAM
C
      SUBROUTINE OUTNUM(R)
C
      DOUBLE PRECISION R
      INTEGER VAL(40),PTR
      PTR=1
      CALL HPXA1(VAL,38,PTR,9,R,I)
      IF (PTR.EQ.0) PTR=41
      VAL(38)='.'
      VAL(39)='.'
      VAL(40)='.'
      PTR=PTR-1
      WRITE (3,1) (VAL(I),I=1,PTR)
1     FORMAT (30X,40A1)
      RETURN
C
      END
C HPXA1 - MINCED SEMXA1 FOR USE WITH HP CALCULATOR
C
      SUBROUTINE HPXA1 (BUFFER,LENGTH,PTR,FUNCT,VALUE,IVALUE)
C
      DOUBLE PRECISION VALUE,V,SIGN,THRESH,BASE,Pi
      LOGICAL EXPON,FRPT,NST,A1FORM,SEMLU
      INTEGER BUFFER(1),PTR,S(4),FUNCT,CH,EX,COL,S1,S2,S3
      INTEGER SPACE,POINT,E,PLUS,A1(64),SHARP
      EQUIVALENCE (S(1),S1),(S(2),S2),(S(3),S3)
C
      COMMON /HPCALC/ Pi,IBASE,IDIGIT,STACK(64)
      DATA SPACE/32/,POINT/46/,E/5/,PLUS/43/,MINUS/45/,SHARP/35/
      DATA S(4)/0/
C
C THIS CHARACTER TABLE SHOULD BE ADJUSTED TO REFLECT ANY LOCAL
C CHARACTER REASSIGNMENTS
      DATA A1/1H@,1HA,1HB,1HC,1HD,1HE,1HF,1HG,
     1   1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,
     2   1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,
     3   1HX,1HY,1HZ,1H[,1H\,1H],1H^,1H_,
     4   1H ,1H!,1H",1H#,1H$,1H%,1H&,1H',
     5   1H(,1H),1H*,1H+,1H,,1H-,1H.,1H/,
     6   1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,
     7   1H8,1H9,1HA,1HB,1HC,1HD,1HE,1HF/
C
      A1FORM=FUNCT.GE.8
      BASE=IBASE
C
C FUNCTIONS 4/9: WRITE NUMBER
  260 V=VALUE
  270 EX=0
      IF (V) 271,320,280
  271 V=-V
      I=MINUS
      IGO=2
      GOTO 1
  280 SIGN=BASE*BASE*BASE*BASE*BASE*BASE*BASE
  290 IF (V.LT.SIGN) GOTO 300
      V=V/BASE
      SIGN=BASE
      EX=EX+1
      GOTO 290
  300 SIGN=1./BASE/BASE/BASE
  310 IF (V.GE.SIGN) GOTO 320
      V=V*BASE
      SIGN=1.
      EX=EX-1
      GOTO 310
C ABS(VALUE)=V*BASE**EX
  320 THRESH=V*1E-15
      V=V+THRESH*.5
      SIGN=BASE*BASE*BASE*BASE*BASE*BASE
      COL=0
      NST=.TRUE.
  330 COL=COL+1
      K=V/SIGN
      IF (COL.GT.8) GOTO 380
      IF (COL-7) 350,360,340
  340 IF (V.LE.THRESH) GOTO 390
      I=POINT
      IGO=3
      GOTO 1
  350 IF (K.EQ.0.AND.NST) GOTO 370
360   I=K+48
      IGO=4
      GOTO 1
365   NST=.FALSE.
      V=V-K*SIGN
  370 SIGN=SIGN/BASE
      GOTO 330
  380 IF (V.GT.THRESH) GOTO 360
  390 IF (EX.EQ.0) RETURN
      I=E
      IGO=5
      GOTO 1
395   V=EX
      GOTO 270
  440 FORMAT (30X,200A1)
C
C
C BUFFER DEPOSIT ROUTINE
1     IF (A1FORM) I=A1(I+1)
      BUFFER(PTR)=I
      PTR=PTR+1
      IF (PTR.GT.LENGTH) GOTO 40
      GOTO (280,280,350,365,395), IGO
40    PTR=0
      RETURN
C
C COPYRIGHT:  IMAGE TECHNIQUES OF CAMBRIDGE LTD 1978
C
      END

      Double precision function ALIMIT(Dr0,Dr)
C
C  Subroutine to compute resolution limit (in A) due to effect of
C  incident beam convergence.
C  Ref: O'Keefe & Pitt, Electron Microscopy 1980, vol.1, pp.122-123.
C
      Double precision Pi,STACK,Dr0,Dr(3),Afact,Bfact,Cfact,Dfact,Efact
      Double precision Alpha,Wavl,Cs,Defocus,Splus,Sminus,Dzero
      COMMON /HPCALC/ Pi,IBASE,IDIGIT,STACK(64)
C
      Dzero=0.0D0
C
C  Assume Alpha (milliradians) is on top on stack
C  thus Alpha MUST be the last entry --
      Alpha = Dr0 * 1.0D-3
C
C  Next three on stack are Wavl, Cs and Defocus
C  (not nec. in that order)
C
Case where Defocus is negative or zero
C  If any one of Dr(i) is negative (or zero),
C  then it MUST be Defocus --
      i = 0
10    i = i+1
      If (i.lt.4) then
        If (Dr(i).le.Dzero) then
          Defocus = Dr(i)
          Do 20 j=i,2
            Dr(j) = Dr(j+1)
20        continue
          Wavl = Dmin1(Dr(1),Dr(2))
          Cs = 1.0D7*Dmax1(Dr(1),Dr(2))
          Goto 100
        Else
          goto 10
        EndIf
      EndIf
C
Case where Defocus is NOT negative or zero --
C  But Defocus is probably larger than Wavl and Cs(mm)
C  whereas Wavl will be smaller than Defocus and Cs(mm).
C  So test for largest and smallest --
      iDef = 1
      iWav = 2
      iCs = 3
      Wavl = Dr(1)
      Defocus = Dr(2)
      Do 30 i=1,3
        If (Dr(i).ge.Defocus) then
          Defocus = Dr(i)
          iDef = i
        Else if (Dr(i).le.Wavl) then
          Wavl = Dr(i)
          iWav = i
        EndIf
30    Continue
      Do 40 i=1,3
        if (i.ne.iDef.and.i.ne.iWav) iCs=i
40    Continue
C      write (3,1200) iWav,iDef,iCs
C1200  format('0iWav,iDef,iCs =',3i5)
      Cs = 1.0D07*Dr(iCs)
C
C  Now we have Alpha,Wavl,Cs,Defocus
100   Continue
C      write(3,1000) Alpha,Wavl,Cs,Defocus
C1000  format(' Alpha,Wavl,Cs,Defocus =',/,x,4f16.5)
C
C  Construct factors
      Afact = Defocus/3.0D0
      Bfact = 3.3D0/(4.0D0*pi*Alpha)
      Cfact = wavl*wavl*Cs
      Dfact = Afact*Afact*Afact/Cfact + Bfact*Bfact
      If (Dfact.le.D0) then
        Goto 9000
      Else
        Efact = Dsqrt(Dfact)
      EndIf
      Splus = (Bfact+Efact)/Cfact
      Sminus = (Bfact-Efact)/Cfact
      Alimit = Dcubert(Splus) + Dcubert(Sminus)
      Alimit = 1.0D0/Alimit
C      write(3,1100) Afact,Bfact,Cfact,Dfact,Efact,Splus,Sminus,Alimit
C1100  format(' A,B,C,D,S+,S-,LIMIT=', 4F15.4,/,21x,4f15.4)
      Return
9000  Alimit = D0
      Return
      END

      Double precision function Dcubert(Dx)
C  To compute CUBE ROOT in double precision
      Double precision Dx,DAx,D0,Dsmall
      Data D0,Dsmall/0.0D0,1.0D-18/
      Dcubert = D0
      DAx = DABS(Dx)
      If (DAx.gt.Dsmall) then
        Dcubert = DEXP(DLOG(DABS(Dx))/3)
        If (Dx.lt.D0) then
          Dcubert=-Dcubert
        EndIf
      EndIf
      Return
      End