C   ******************** SEMITIP_V3 ************************
C
C   CALLING PROGRAM FOR E-FIELD COMPUTATIONS FOR TIP-SEMICONDUCTOR
C   TAKES INPUT FROM FILE FORT.9
C
      DIMENSION VAC(2,2048,2048),SEM(2,2048,2048),VSINT(2,2048),
     &R(2048),S(2048),DELV(2048),ITMAX(10),EP(10),RHOC(50000),
     &RHOS(50000)
      LOGICAL TIP (2048,2048)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV
      COMMON/PROTRU/RAD2
      COMMON/SURF/EN,DENS
      DATA EPSIL0/8.854185E-12/E/1.60210E-19/
      PI=4.*ATAN(1.)
C
      NRDIM=2048
      NVDIM=2048
      NSDIM=2048
      NE=50000
C
C   SEMICONDUCTOR AND TIP PARAMETERS
C
      READ(9,*) NPARM
      DO 900 IPARM=1,NPARM
      READ(9,*) SLOPE
      READ(9,*) SEP
      READ(9,*) RAD
      READ(9,*) RAD2
      READ(9,*) BIAS
      READ(9,*) DELPHI
      PHITIP=BIAS+DELPHI
      READ(9,*) CD
      READ(9,*) CA
      THETA=360.*ATAN(1./SLOPE)/PI
      WRITE(6,*) ' '
      WRITE(16,*) ' '
      WRITE(6,*) 'SEP, BIAS =',SEP,BIAS
      WRITE(16,*) 'SEP, BIAS =',SEP,BIAS
      WRITE(6,*) 'CONTACT POTENTIAL, TIP POTENTIAL =',DELPHI,PHITIP
      WRITE(16,*) 'CONTACT POTENTIAL, TIP POTENTIAL =',DELPHI,PHITIP
      WRITE(6,*) 'DOPING =',CD,CA
      WRITE(16,*) 'DOPING =',CD,CA
      WRITE(6,*) 'RAD, SLOPE, ANGLE =',RAD,SLOPE,THETA 
      WRITE(16,*) 'RAD, SLOPE, ANGLE =',RAD,SLOPE,THETA 
      READ(9,*) EGAP
      READ(9,*) ED
      READ(9,*) EA
      READ(9,*) ACB
      READ(9,*) AVB
      READ(9,*) EPSIL
      READ(9,*) TEM
      TK=TEM*8.617E-5
      READ(9,*) IDEG
      READ(9,*) IINV
      READ(9,*) DENS
      READ(9,*) EN
      WRITE(6,*) 'SURFACE STATE DENSITY, EN =',DENS,EN
      WRITE(16,*) 'SURFACE STATE DENSITY, EN =',DENS,EN
C
C   FIND FERMI-LEVEL POSITION
C
      CALL EFFIND(EF)
      WRITE(6,*) 'DOPING, FERMI-LEVEL =',CD,CA,EF
      WRITE(16,*) 'DOPING, FERMI-LEVEL =',CD,CA,EF
C
C   CHARGE NEUTRALITY LEVEL = EN
C
      PHIS=EF-EN
C
C   1-D SOLUTION FOR BAND BENDING
C
      IF ((CD-CA).EQ.0.) GO TO 105
      W=1.E9*SQRT(2.*EPSIL*EPSIL0*AMAX1(1.,ABS(PHITIP))/
     &   (ABS(CD-CA)*1.E6*E))
      WRITE(6,*) '1-D ESTIMATE OF DEPLETION WIDTH (NM) =',W
      WRITE(16,*) '1-D ESTIMATE OF DEPLETION WIDTH (NM) =',W
      GO TO 110
105   W=1.E10
C
C   CONSTRUCT TABLES OF CHARGE DENSITY VALUES
C
110   ESTART=AMIN1(EF,EF-PHITIP,EF-PHIS)
      EEND=AMAX1(EF,EF-PHITIP,EF-PHIS)
      ETMP=EEND-ESTART
      ESTART=ESTART-2.*ETMP
      EEND=EEND+2.*ETMP
      DELE=(EEND-ESTART)/FLOAT(NE-1)
      WRITE(6,*) 'ESTART,EEND,NE =',ESTART,EEND,NE
      WRITE(16,*) 'ESTART,EEND,NE =',ESTART,EEND,NE
      CALL SEMIRHO(DELE,ESTART,NE,RHOC)
      CALL SURFRHO(DELE,ESTART,NE,RHOS)
C
C   SEMICONDUCTOR GRID SIZE AND SPACING
C
      READ(9,*) NR
      READ(9,*) NV
      READ(9,*) NS
      READ(9,*) SIZE
      DELR=RAD*SIZE
      IF (RAD2.NE.0.) DELR=AMIN1(RAD2,DELR)
      DELR=AMIN1(DELR,W/NR)
C     NV=MAX0(NV,NINT(RAD*SLOPE**2/SEP))
      DELS=RAD*SIZE
      IF (RAD2.NE.0.) DELS=AMIN1(RAD2,DELS)
      DELS=AMIN1(DELS,W/NS)
C
C   ITERATION LIMITS
C
      READ(9,*) IPMAX
      READ(9,*) (ITMAX(IP),IP=1,IPMAX)
      READ(9,*) (EP(IP),IP=1,IPMAX)
      READ(9,*) IWRIT
C
C   SOLVE THE PROBLEM !
C
      CALL SEMITIP2(SEP,RAD,SLOPE,ETAT,A,Z0,C,VAC,TIP,SEM,VSINT,
     &R,S,DELV,DELR,DELS,NRDIM,NVDIM,NSDIM,NR,NV,NS,PHITIP,IWRIT,
     &NE,ESTART,DELE,RHOC,RHOS,EF,ITMAX,EP,IPMAX,PHI0,IERR)
      WRITE(6,500) NR,NS,NV,IERR
      WRITE(16,500) NR,NS,NV,IERR
500   FORMAT(' RETURN FROM SEMTIP2, NR,NS,NV,IERR =',4I7)
C
C   OUTPUT
C
      WRITE(10,510) RAD,SEP,BIAS,THETA,PHI0
510   FORMAT(' ',5G12.4)
C
C   PLOT CROSS-SECTIONAL PROFILE
C
      IF (IWRIT.GE.1) THEN
      DO 705 J=NV,1,-1
        WRITE(11,*) -J*DELV(1),VAC(1,1,J)
705   CONTINUE
      WRITE(11,*) 0.,VSINT(1,1)
      DO 710 J=1,NS-1
          WRITE(11,*) S(J),SEM(1,1,J)
710   CONTINUE
C
C   PLOT SURFACE POTENTIAL
C
      WRITE(12,*) -R(1),VSINT(1,1)
      DO 720 I=1,NR-1
          WRITE(12,*) R(I),VSINT(1,I)
720   CONTINUE
      END IF
C
C   PLOT CONTOURS
C
      IF (IWRIT.GE.2) THEN
      READ(9,*) NUMC
      READ(9,*) DELP
      CALL CONTR(ETAT,VAC,TIP,SEM,VSINT,R,S,DELV,NRDIM,NVDIM,
     &NSDIM,NR,NV,NS,NUMC,DELP)
      END IF
C
C   OUTPUT ENTIRE POTENTIAL
C
      IF (IWRIT.GE.3) THEN
      nrecl=36+4*nr*(nv+ns+3)+4*ns
      open(unit=13,file='fort.13',access='direct',recl=nrecl)                                            
      write(13,rec=1) nr,nv,ns,sep,rad,rad2,slope,bias,epsil,
     &((vac(1,i,j),i=1,nr),j=1,nv),
     &((sem(1,i,j),i=1,nr),j=1,ns),(vsint(1,i),i=1,nr),
     &(r(i),i=1,nr),(s(j),j=1,ns),(delv(i),i=1,nr)
      END IF
C
 900  CONTINUE
      WRITE(6,*) 'PRESS THE ENTER KEY TO EXIT'
      WRITE(16,*) 'PRESS THE ENTER KEY TO EXIT'
      READ(5,*)
C
      stop
      end
C
      real function p(r)
      common/protru/rad2
      p=0.
      if (r.lt.rad2) p=sqrt(rad2**2-r**2)
      return
      end
c
      real function sig(EF,PHI)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      common/SURF/EN,DENS
      PI=4.*ATAN(1.)
c
c     width=fwhm/(2.*sqrt(2.*alog(2.)))
c     sig=-exp(-(EF-PHI-(EN+ecent))**2/(2.*width**2))+
c     &     exp(-(EF-PHI-(EN-ecent))**2/(2.*width**2))
c     sig=sig*dens/(SQRT(2.*pi)*width)
c     sig=exp(-(EF-PHI-EN)**2/(2.*width**2))
c     sig=sig*dens/(SQRT(2.*pi)*width)
c
c     sig=0.
c     if ((ef-phi).lt.0.) return
c     if ((ef-phi).gt.(eginv-eaa)) return
c     sig=exp(-(EF-PHI)/width)-
c     &    exp((EF-PHI-(eginv-eaa))/width)
c     sig=sig*dens/WIDTH
c
      sig=dens
      if ((ef-phi).gt.en) sig=-sig
      return
      end