C   ******************** POTINT ************************
C
C   INTEGRATE CHARGE DENSITIES TO GET TIP-INDUCED BAND BENDING
C   CURRENT PROGRAM WORKS FOR ZERO DOPING AND IDEAL TIP
C
C   VERSION 6.0 - MAR/11, DERIVED FROM VERY OLD VERSION OF PROGRAM   
C
      FUNCTION POTINT(SEP,ETAT,VAC,TIP,SEM,VSINT,R,DELR,DELV,DELS,
     &NR,NV,NS,NRDIM,NVDIM,NSDIM,BIAS,EPSIL)
C
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &R(NRDIM),DELR(NRDIM),DELV(NRDIM),DELS(NSDIM)
      LOGICAL TIP(NRDIM,NVDIM)
      DOUBLE PRECISION SUM
C
      SUM=0.D0
      TERM=0.D0
      RR=0.D0
      DO 400 I=1,NR-1
         DO 210 J=NV,1,-1
            IF (.NOT.TIP(I,J)) GO TO 220
210      CONTINUE
220      JSAV=J
         DPHI=BIAS-VAC(2,I,JSAV)
         RRSAV=RR
         RR=R(I)*(SQRT(1.-(ETAT*(JSAV+1)/NV)**2)+
     &            SQRT(1.-(ETAT*JSAV/NV)**2))/2.
         DR=R(I)*(SQRT(1.-(ETAT*(JSAV+1)/NV)**2)-
     &            SQRT(1.-(ETAT*JSAV/NV)**2))
         ZZ=(JSAV+0.5)*DELV(I)
         DZ=DELV(I)
         SIG=DPHI/SQRT(DR**2+DZ**2)
         B=RR*ETAT**2/(ZZ*(1.-ETAT**2))
         TERMSAV=TERM
C   FIRST TERM BELOW IS POTENTIAL AT SURFACE RELATIVE TO TIP APEX;
C   SECOND TERM IS POTENTIAL AT DISTANCE SEP INSIDE TIP RELATIVE TO TIP APEX
         TERM=RR*SIG*(SQRT(1.+B**2))*0.5*(1./SQRT(RR**2+(ZZ-SEP)**2)-
     &        (EPSIL-1.)/(SQRT(RR**2+(ZZ+SEP)**2)*(EPSIL+1.))-
     &        2./(SQRT(RR**2+ZZ**2)*(EPSIL+1.)))
C        TERM=RR*SIG*(SQRT(1.+B**2))*0.5*(1./SQRT(RR**2+(ZZ-2.*SEP)**2)-
C     &        (EPSIL-1.)/(SQRT(RR**2+(ZZ+2.*SEP)**2)*(EPSIL+1.))-
C     &        1./SQRT(RR**2+(ZZ-2.*SEP)**2)+
C     &        (EPSIL-1.)/(SQRT(RR**2+(ZZ+2.*SEP)**2)*(EPSIL+1.)))
         SUM=SUM+(TERM+TERMSAV)*(RR-RRSAV)/2.
400   CONTINUE
      POTINT=SUM
      RETURN
      END