C   ******************** SEMITIP2 ************************
C
C   E-FIELD COMPUTATIONS FOR TIP IN PROXIMITY TO SEMICONDUCTOR,
C   WITH CIRCULAR SYMMETRY
C   VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, DEC 2002
C
C   SEMI COMMON BLOCK:
C      EGAP=BAND GAP (EV)
C      ED=DONOR BINDING ENERGY (EV)
C      EA=ACCEPTOR BINDING ENERGY (EV)
C      ACB=CONDUCTION BAND EFFECTIVE MASS
C      AVB=VALENCE BAND EFFECTIVE MASS
C      CD=DONOR DENSITY (CM^-3)
C      CA=ACCEPTOR DENSITY (CM^-3)
C      EPSIL=DIELECTRIC CONSTANT OF SEMICONDUCTOR
C      TK=TEMPERATURE (EV)
C      IDEG=1 IF DEGENERATE DOPING, 0 OTHERWISE
C
C   INPUT AND OUTPUT PARAMETERS:
C      SEP=SEPARATION BETWEEN SEMICONDUCTOR AND END OF TIP (NM)
C      RAD=RADIUS OF HYPERBOLOID WHICH FORMS TIP (NM)
C      SLOPE=SLOPE OF TIP SHANK
C      ETAT=ETA PARAMETER FOR THE HYPERBOLIC TIP
C      AT=LOCATION OF FOCUS FOR THE HYPERBOLIC TIP
C      ETA1=ETA PARAMETER FOR THE HYPERBOLIC COORDINATE SYSTEM
C      A1=A2=LOCATION OF FOCUS FOR THE HYPERBOLIC COORDINATE SYSTEM
C      VAC=ARRAY OF POTENTIAL VALUES IN VACUUM
C      SEM=ARRAY OF POTENTIAL VALUES IN SEMICONDUCTOR
C      VSINT=ARRAY OF POTENTIAL VALUES AT SEMICONDUCTOR SURFACE
C      DELV=ARRAY OF Z-SPACING VALUES IN VACUUM
C      DELR=R-SPACING IN VACUUM AND SEMICONDUCTOR
C      DELS=Z-SPACING IN SEMICONDUCTOR
C      NRDIM=R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
C      NVDIM=Z-DIMENSION FOR VAC ARRAY
C      NSDIM=Z-DIMENSION FOR SEM ARRAY
C      NR=CURRENT VALUE OF R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
C      NV=CURRENT VALUE OF Z-DIMENSION FOR VAC ARRAY
C      NS=CURRENT VALUE OF Z-DIMENSION FOR SEM ARRAY
C      BIAS=BIAS OF TIP RELATIVE TO A POINT FAR INSIDE THE SEMICONDUCTOR
C      IWRIT=1 IF WANT WRITTEN OUTPUT, 0 OTHERWISE
C      NE=NUMBER OF POINTS IN CHARGE DENSITY TABLE
C      ESTART=STARTING ENERGY IN TABLE
C      DELE=ENERGY STEP SIZE IN TABLE
C      RHOC=CHARGE DENSITY TABLE
C      EF=FERMI-LEVEL POSITION
C      ITMAX=ARRAY OF MAXIMUM ITERATION SPECIFICATIONS
C      EP=ARRAY OF STOPPING CRITERIONS
C      IPMAX=MAXIMUM NUMBER OF SCALING LOOPS FOR GRID DOUBLING
C      PHI0=MAX SURFACE BAND BENDING
C      IERR=ERROR INDICATION (=1 IF GRID SIZE TOO SMALL, =2 IF UNSTABLE SOLN)
C
C
C   SEMMIN USED TO FIND OPTIMAL UPDATE EQUATION FOR SOLN IN SEMICONDUCTOR
C
      FUNCTION SEMMIN(PHI,RHOC,NE,IERR)
      DIMENSION RHOC(NE)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/STIP/BIAS,DELR,DELS,DELETA,A1,ESTART,DELE,EF,STEMP
      DATA EEP/1.80943E-20/
      ENER=EF-PHI
      IENER=NINT((ENER-ESTART)/DELE)+1
      IF (IENER.GE.1.AND.IENER.LE.NE) GO TO 100
         WRITE(6,*) 'ERROR - ENERGY INDEX OUT OF BOUNDS'
         WRITE(16,*) 'ERROR - ENERGY INDEX OUT OF BOUNDS'
         IERR=2
         IF (IENER.LT.1) IENER=1
         IF (IENER.GT.NE) IENER=NE
100   TEMP=STEMP-RHOC(IENER)*EEP/EPSIL
      SEMMIN=ABS(PHI-TEMP/(2./DELR**2+2./DELS**2))
      RETURN
      END
C
C   GOLDEN SECTION SEARCH ROUTINE (USED WITH SEMMIN)
C
      SUBROUTINE GSECT2(XMIN,XMAX,EP,RHOC,NE,IERR)
      DIMENSION RHOC(NE)
      DATA GS/0.3819660/
      IF (XMAX.EQ.XMIN) RETURN
      IF (EP.EQ.0.) RETURN
      IF (XMAX.LT.XMIN) THEN
         TEMP=XMAX
         XMAX=XMIN
         XMIN=TEMP
         END IF
      DELX=XMAX-XMIN
      XA=XMIN+DELX*GS
      FA=SEMMIN(XA,RHOC,NE,IERR)
      XB=XMAX-DELX*GS
      FB=SEMMIN(XB,RHOC,NE,IERR)
 100  DELXSAV=DELX
      IF (DELX.LT.EP) RETURN
      IF (FB.LT.FA) GO TO 200
      XMAX=XB
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XB=XA
      FB=FA
      XA=XMIN+DELX*GS
      FA=SEMMIN(XA,RHOC,NE,IERR)
      GO TO 100
 200  XMIN=XA
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XA=XB
      FA=FB
      XB=XMAX-DELX*GS
      FB=SEMMIN(XB,RHOC,NE,IERR)
      GO TO 100
      END
C
C   MAIN PROGRAM
C
      SUBROUTINE SEMITIP2(SEP,RAD,SLOPE,ETAT,AT,ETA1,A2,VAC,SEM,VSINT,
     &DELV,DELR1,DELS1,NRDIM,NVDIM,NSDIM,NR1,NV1,NS1,BIAS1,IWRIT,
     &NE,ESTART1,DELE1,RHOC,EF1,ITMAX,EP,IPMAX,PHI0,IERR)
C
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &DELV(NRDIM),ITMAX(10),EP(10),RHOC(NE)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/SIZE/NR,NV,NS
      COMMON/STIP/BIAS,DELR,DELS,DELETA,A1,ESTART,DELE,EF,STEMP
      BIAS=BIAS1
      ESTART=ESTART1
      DELE=DELE1
	EF=EF1
      NR=NR1
      NV=NV1
      NS=NS1
	DELR=DELR1
	DELS=DELS1
      IERR=0
C
C   CONSTRUCT THE TIP AND VACUUM GRID AND INITIALIZE POTENTIAL
C
      ETAT=1./SQRT(1.+1./SLOPE**2)
      AT=RAD*SLOPE**2/ETAT
      IF (SEP.GT.AT*ETAT) GO TO 50
C   CASE 1 - BASIS TIP HAS SAME ETA & A AS ACTUAL TIP, BUT DIFFERENT SEPARATION
         ETA1=ETAT
         A1=AT
         A2=A1
         GO TO 60
C   CASE 2 - BASIS TIP HAS SAME SEPARATION AS ACTUAL TIP, BUT DIFFERENT ETA & A
50       ETA1=1./SQRT(1.+RAD/SEP)
         A1=SQRT(SEP*(RAD+SEP))
60      CONTINUE
      DELETA=ETA1/FLOAT(NV)
      CETA1=ALOG((1.+ETA1)/(1.-ETA1))
      IF (IWRIT.EQ.0) GO TO 100
         WRITE(6,*) 'ACTUAL ETA, A, S =',ETAT,AT,SEP
         WRITE(16,*) 'ACTUAL ETA, A, S =',ETAT,AT,SEP
         WRITE(6,*) 'COORD  ETA, A, S =',ETA1,A1,ETA1*A1
         WRITE(16,*) 'COORD  ETA, A, S =',ETA1,A1,ETA1*A1
         WRITE(6,*) 'NR,NS,NV =',NR,NS,NV
         WRITE(16,*) 'NR,NS,NV =',NR,NS,NV
         WRITE(6,*) 'DELR,DELS,DELV =',DELR,DELS,A1*DELETA
         WRITE(16,*) 'DELR,DELS,DELV =',DELR,DELS,A1*DELETA
100   CONTINUE
      IF (NR.GT.NRDIM.OR.NV.GT.NVDIM.OR.NS.GT.NSDIM) GO TO 600
      DO 150 I=1,NR
         X0=DELR*(I-0.5)
         DELV(I)=SQRT(A1**2+X0**2)*ETA1/FLOAT(NV)
         DO 120 J=1,NV-1
            XSI=SQRT(1.+(X0/A1)**2)
            ETA=J*ETA1/FLOAT(NV)
            Z=A1*XSI*ETA
            R=A1*SQRT((XSI**2-1.)*(1.-ETA**2))
            IF (Z.GT.(ETAT*(SQRT(R**2/(1.-ETAT**2)+AT**2)-AT)+SEP-P(R)))
     &         GO TO 110
            VAC(1,I,J)=0.
            VAC(2,I,J)=0.
            GO TO 120
110         VAC(1,I,J)=BIAS
            VAC(2,I,J)=BIAS
120      CONTINUE
         VAC(1,I,NV)=BIAS
         VAC(2,I,NV)=BIAS
150   CONTINUE
      DO 160 I=1,NR
         VSINT(1,I)=0.
         VSINT(2,I)=0.
160   CONTINUE
      DO 180 J=1,NS
         DO 170 I=1,NR
            SEM(1,I,J)=0.
            SEM(2,I,J)=0.
170      CONTINUE
180   CONTINUE
      IF (DELV(1).LT.SEP) GO TO 190
         IERR=1
         WRITE(6,*) 'ERROR - VACUUM GRID SPACING TOO LARGE'
         WRITE(16,*) 'ERROR - VACUUM GRID SPACING TOO LARGE'
190   CONTINUE
      IF (IERR.NE.0) GO TO 600
C
C   INITIAL GUESS
C
      DO 250 I=1,NR
         DO 210 J=NV,1,-1
            IF (VAC(1,I,J).NE.BIAS) GO TO 220
210      CONTINUE
220      JSAV=J
         DO 230 J=JSAV,1,-1
            ETA=J*DELETA*FLOAT(NV)/FLOAT(JSAV+1)
            VAC(1,I,J)=BIAS*ALOG((1.+ETA)/(1.-ETA))/CETA1
            VAC(2,I,J)=VAC(1,I,J)
230      CONTINUE
250   CONTINUE
      PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
C
C   ITERATE POISSON'S EQN
C
      DO 590 IP=1,IPMAX
C
         IF (IWRIT.NE.0) THEN
	      WRITE(6,*) 'SOLUTION #',IP
	      WRITE(16,*) 'SOLUTION #',IP
	   END IF
         ITM=ITMAX(IP)
         EPI=EP(IP)
         CALL ITER2(VAC,SEM,VSINT,DELV,NRDIM,NVDIM,NSDIM,
     &              EPI,ITM,PHI0,IWRIT,RHOC,NE,IERR)
         IF (IP.EQ.IPMAX) GO TO 600
C
         IF (NR*2.GT.NRDIM.OR.NV*2.GT.NVDIM.OR.NS*2.GT.NSDIM.OR
     &      .IERR.NE.0) GO TO 500
C
C   DOUBLE ALL ARRAYS
C
         NR=NR*2
         NS=NS*2
         NV=NV*2
         DELR=DELR/2.
         DELS=DELS/2.
         DELETA=DELETA/2.
         IF (IWRIT.EQ.0) GO TO 300
            WRITE(6,*) 'NR,NS,NV =',NR,NS,NV
            WRITE(16,*) 'NR,NS,NV =',NR,NS,NV
            WRITE(6,*) 'DELR,DELS,DELV =',DELR,DELS,A1*DELETA
            WRITE(16,*) 'DELR,DELS,DELV =',DELR,DELS,A1*DELETA
300      CONTINUE
         DO 450 I=1,NR
            X0=DELR*(I-0.5)
            DELV(I)=SQRT(A1**2+X0**2)*ETA1/FLOAT(NV)
            DO 430 J=1,NV-1
               X2M1=(X0/A1)**2
               XSI=DSQRT(1.D0+X2M1)
               ETA=J*ETA1/FLOAT(NV)
               Z=A1*XSI*ETA
               X=A1*SQRT(X2M1*(1.-ETA**2))
               IF (Z.GT.(ETAT*(SQRT(X**2/(1.-ETAT**2)+AT**2)-AT)+SEP
     &            -P(X))) GO TO 410
               IF (J.EQ.1) GO TO 390
               IF (MOD(J,2).NE.0) GO TO 380
               TEMP=VAC(2,(I-1)/2+1,J/2)
               GO TO 400
380            TEMP=(VAC(2,(I-1)/2+1,J/2)+VAC(2,(I-1)/2+1,(J+1)/2))/2.
               GO TO 400
390            TEMP=(VAC(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2.
400            IF (TEMP.EQ.BIAS) TEMP=0.9999*TEMP
               VAC(1,I,J)=TEMP
               GO TO 430
410         VAC(1,I,J)=BIAS
               VAC(2,I,J)=BIAS
430         CONTINUE
            VAC(1,I,NV)=BIAS
450      CONTINUE
         DO 460 I=1,NR
            VSINT(1,I)=VSINT(2,(I-1)/2+1)
460      CONTINUE
         DO 465 I=1,NR
             SEM(1,I,1)=(SEM(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2.
465      CONTINUE
         DO 475 J=2,NS,2
            DO 470 I=1,NR
                SEM(1,I,J)=SEM(2,(I-1)/2+1,J/2)
470         CONTINUE
475      CONTINUE
         DO 485 J=3,NS,2
            DO 480 I=1,NR
                SEM(1,I,J)=(SEM(2,(I-1)/2+1,J/2)+
     &                     +SEM(2,(I-1)/2+1,(J+1)/2))/2.
480         CONTINUE
485      CONTINUE
         DO 492 J=1,NV
            DO 490 I=1,NR
               VAC(2,I,J)=VAC(1,I,J)
490         CONTINUE
492      CONTINUE
         DO 494 I=1,NR
            VSINT(2,I)=VSINT(1,I)
494      CONTINUE
         DO 498 J=1,NS
            DO 496 I=1,NR
               SEM(2,I,J)=SEM(1,I,J)
496         CONTINUE
498      CONTINUE
         GO TO 590
C
500      IF (NS*2.GT.NSDIM.OR.IERR.NE.0) GO TO 600
C
C   DOUBLE ONLY SEM ARRAY ELEMENTS
C
         NS=NS*2
         DELS=DELS/2.
         IF (IWRIT.EQ.0) GO TO 505
            WRITE(6,*) 'NS =',NS
            WRITE(16,*) 'NS =',NS
            WRITE(6,*) 'DELS =',DELS
            WRITE(16,*) 'DELS =',DELS
505      CONTINUE
         DO 510 I=1,NR
             SEM(1,I,1)=(SEM(2,I,1)+VSINT(2,I))/2.
510      CONTINUE
         DO 520 J=2,NS,2
            DO 515 I=1,NR
                SEM(1,I,J)=SEM(2,I,J/2)
515         CONTINUE
520      CONTINUE
         DO 530 J=3,NS,2
            DO 525 I=1,NR
                SEM(1,I,J)=(SEM(2,I,J/2)+
     &                     +SEM(2,I,(J+1)/2))/2.
525         CONTINUE
530      CONTINUE
         DO 550 J=1,NS
            DO 540 I=1,NR
               SEM(2,I,J)=SEM(1,I,J)
540         CONTINUE
550      CONTINUE
590   CONTINUE
C
600   CONTINUE
      IF (IWRIT.NE.0) THEN
	   WRITE(6,*) 'NUMBER OF ITERATIONS = ',ITM
	   WRITE(16,*) 'NUMBER OF ITERATIONS = ',ITM
         WRITE(6,*) 'BAND BENDING AT MIDPOINT = ',PHI0
         WRITE(16,*) 'BAND BENDING AT MIDPOINT = ',PHI0
	END IF
C
      NR1=NR
      NV1=NV
      NS1=NS
      DELR1=DELR
      DELS1=DELS
      RETURN
      END
C
      INTEGER FUNCTION IADS(I)
      COMMON/SIZE/NR,NV,NS
      IF (I.EQ.0) GO TO 10
      IF (I.GT.NR) GO TO 20
      IADS=I
      RETURN
10    IADS=1
      RETURN
20    IADS=NR
      RETURN
      END
C
C   3-D SOLUTION FOR POTENTIAL
C
      SUBROUTINE ITER2(VAC,SEM,VSINT,DELV,NRDIM,NVDIM,NSDIM,EP,ITMAX,
     &PHI0,IWRIT,RHOC,NE,IERR)
C
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &DELV(NRDIM),RHOC(NE)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,DEGEN
      COMMON/SIZE/NR,NV,NS
      COMMON/STIP/BIAS,DELR,DELS,DELETA,A1,ESTART,DELE,EF,STEMP
      DATA EEP/1.80943E-20/
C
      DO 500 ITER=1,ITMAX
         DO 360 I=1,NR
            X2M1=((I-0.5)*DELR/A1)**2
            XSI=DSQRT(1.D0+X2M1)
            IF (X2M1.LT.1.E-4) GO TO 320 
                DELXSI=DSQRT(1.D0+(I*DELR/A1)**2)-
     &                DSQRT(1.D0+((I-1)*DELR/A1)**2)
               GO TO 330
320            DELXSI=(I-0.5)*(DELR/A1)**2
330         CONTINUE
            DO 350 J=1,NV-1
               IF (VAC(1,I,J).EQ.BIAS) GO TO 360
               ETA=J*DELETA
               IF (J.EQ.1) GO TO 335
                  TEMP=X2M1*(VAC(1,IADS(I+1),J)+VAC(1,IADS(I-1),J))
     &               /DELXSI**2+
     &              (1.-ETA**2)*(VAC(1,I,J+1)+VAC(1,I,J-1))/DELETA**2+
     &              (VAC(1,IADS(I+1),J)-VAC(1,IADS(I-1),J))*XSI/DELXSI-
     &              (VAC(1,I,J+1)-VAC(1,I,J-1))*ETA/DELETA
               GO TO 340
335               TEMP=X2M1*(VAC(1,IADS(I+1),1)+VAC(1,IADS(I-1),1))
     &            /DELXSI**2+
     &           (1.-ETA**2)*(VAC(1,I,2)+VSINT(1,I))/DELETA**2+
     &           (VAC(1,IADS(I+1),1)-VAC(1,IADS(I-1),1))*XSI/DELXSI-
     &           (VAC(1,I,2)-VSINT(1,I))*ETA/DELETA
340            CONTINUE
               VAC(2,I,J)=TEMP/(2.*X2M1/DELXSI**2+
     &                          2.*(1.-ETA**2)/DELETA**2)
350         CONTINUE
360      CONTINUE
         DO 380 J=1,NV
            DO 370 I=1,NR
               VAC(1,I,J)=VAC(2,I,J)
370         CONTINUE
380      CONTINUE
         DO 400 I=1,NR
            VSINT(2,I)=(VAC(1,I,1)/DELV(I)+EPSIL*SEM(1,I,1)/DELS)/
     &                 (1./DELV(I)+EPSIL/DELS)
400      CONTINUE
         DO 410 I=1,NR
            VSINT(1,I)=VSINT(2,I)
410      CONTINUE
         IF (MOD(ITER,100).EQ.0) THEN
            PHISAV=PHI0
            PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
            IF (IWRIT.NE.0) THEN
		     WRITE(6,*) 'ITER,PHI0 =',ITER,PHI0
		     WRITE(16,*) 'ITER,PHI0 =',ITER,PHI0
	      END IF
         END IF
         DO 450 J=1,NS-1
            DO 440 I=1,NR
               R=DELR*SQRT((I-0.5)**2)
               SEMOLD=SEM(1,I,J)
               ENER=EF-SEM(1,I,J)
               IF (J.EQ.1) GO TO 420
                  STEMP=(SEM(1,IADS(I+1),J)+SEM(1,IADS(I-1),J))/DELR**2+
     &             (SEM(1,I,J+1)+SEM(1,I,J-1))/DELS**2+
     &             (SEM(1,IADS(I+1),J)-SEM(1,IADS(I-1),J))/(2.*R*DELR)
	         GO TO 430
420               STEMP=(SEM(1,IADS(I+1),1)+SEM(1,IADS(I-1),1))/DELR**2+
     &             (SEM(1,I,2)+VSINT(1,I))/DELS**2+
     &             (SEM(1,IADS(I+1),1)-SEM(1,IADS(I-1),1))/(2.*R*DELR)
430            CONTINUE
               IENER=NINT((ENER-ESTART)/DELE)+1
               IF (IENER.LT.1.OR.IENER.GT.NE) THEN
                  WRITE(6,*) 'ERROR - ENERGY INDEX OUT OF BOUNDS'
                  WRITE(16,*) 'ERROR - ENERGY INDEX OUT OF BOUNDS'
	         END IF
               IF (IENER.LT.1) IENER=1
               IF (IENER.GT.NE) IENER=NE
               TEMP=STEMP-RHOC(IENER)*EEP/EPSIL
               SEMNEW=TEMP/(2./DELR**2+2./DELS**2)
               NSRCH=10
               DELSEM=(SEMNEW-SEMOLD)/FLOAT(NSRCH)
               CALL GSECT2(SEMOLD,SEMNEW,DELSEM,RHOC,NE,IERR)
               SEM(2,I,J)=(SEMOLD+SEMNEW)/2.
440         CONTINUE
450      CONTINUE
         DO 480 J=1,NS
            DO 470 I=1,NR
               SEM(1,I,J)=SEM(2,I,J)
470         CONTINUE
480      CONTINUE   
         IF ((MOD(ITER,100).EQ.0.AND.ABS(PHI0-PHISAV).LT.EP).OR
     &      .IERR.NE.0) GO TO 510
500   CONTINUE
510   ITMAX=ITER
      PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
      RETURN
      END