      program poissona
c
c   Phi'' = F(Phi,h)
c         = -(e/eps0) sum  Zi ni(Phi,h)
c                barom: ni=  ni(h0) exp(-(mi Phig+Zi e Phi)/k Ti)
c
c      h0 <=  h  <= hN
c
      parameter (ndim=75,ndat=12)
      real*8 Phi(0:ndim),Phicor(0:ndim),hb,h0,hi,hN,Phi0,Phi1,
     &   xi(0:ndim),  D2(ndim,ndim) ,fornb(0:ndim,0:2), PhiNn,Phi0n,
     &   resid(ndim),Phig,one,pi,test,ti,mat(ndim,ndim),densit(ndat),
     &   mass(ndat),temper(ndat),aux,aux1,erfc,den,dden ,fk2
     &  ,fnorm(ndat),PhiN(ndim),PhiNr,denox,hox,charge,aux2,tayl(6)
     &  ,hcrit,err,derr
      real plot(0:500,3),mp,me,scalm(3),scalx(3)
      integer iz(ndat),ic(ndat),ipsi(ndat),ixi(ndat),izeta(ndat)
      character*20 fdata
      character*80 ldata
      character*4 nom(ndat)
      COMMON /cphys/e,ee,Tk,mp,me,RE
      one=1
      pi=4*atan(one)
      spi=sqrt(pi)
      aux=erfc(one/2)
      print *,' test erfc(1/2): ',aux,' err=',one-0.5204998778-aux
      aux=0
      print *,' test 1/2-K2 (0):',fk2(aux)
      aux=1
      print *,' test 1/2-K2 (1):',fk2(aux)
      aux=10
      print *,' test 1/2-K2 (10):',fk2(aux),
     &                exp(-aux*aux)*(aux+one/(2*aux))/spi
      aux=15
      print *,' test 1/2-K2 (15):',fk2(aux),
     &                exp(-aux*aux)*(aux+one/(2*aux))/spi
      print *,'  name of data file ? '
      read *,fdata
      open(unit=1, file=fdata, status='old')
      open(unit=12, file='poissona.ps')
      write(12, '(2A1)')char(37),char(33)
      write(12, '(''/m{moveto}def /l{lineto}def'')')
      write(12, '(''/Courier-Bold findfont 10 scalefont setfont'')')
      read(1,'(a80)')ldata
      read(1,*)izmax,hb,Phi1
      print *,' contents of ',fdata,'  Phi1=',Phi1
      print '(a80)',ldata
      write(12,'(''85 700 m ('',a80,'')show'')')ldata
      izmax=min(izmax,ndat)
      test=0
      denox=0
      hox=7e7
      hcrit=hox
      do i=1,izmax
       read(1,*)nom(i),iz(i),mass(i),densit(i),temper(i),
     &          ic(i),ipsi(i),ixi(i),izeta(i)
       print '(i4,1pe11.4,1x,a4,i4,3e11.4,2x,4i3)',
     &       i,hb,nom(i),iz(i),mass(i),densit(i),temper(i),
     &          ic(i),ipsi(i),ixi(i),izeta(i)
       write(12,'(''85 '',i4,'' m ('',i4,1pe11.4,1x,a4,i4,3e11.4,
     &  2x,4i3,'')show'')') 700-10*i,
     &       i,hb,nom(i),iz(i),mass(i),densit(i),temper(i),
     &          ic(i),ipsi(i),ixi(i),izeta(i)
       test=test+iz(i)*densit(i)
c  SI units & Temp*k:
       if( abs(mass(i)-16).lt. 0.001)denox=densit(i)
       mass(i)=mass(i)*mp
       temper(i)=temper(i)*Tk
      enddo
      print *,' test sum Zi ni = ',test
      write(12,'(''220 '',i4,'' m (test sum Zi ni = '',1p,e12.4,
     &  '' , Phi1= '',0p,f7.2,'')show'')')700-10*(izmax+2),test,Phi1
c alpha, somme des flux et fnorm
      Phig=-9.81*RE/(1+hb/RE)
      sflux=0
      do i=1,izmax
       ali=iz(i)*e*Phi1/( mass(i)*Phig)
c    vinf2
       aux=(-mass(i)*Phig-iz(i)*e*Phi1)/temper(i)
       if(aux.lt.0)aux=-mass(i)*Phig/(temper(i)*1e10)
       aux=sqrt(aux)
       fnorm(i)=(izeta(i)-ipsi(i))*fk2(aux)+ipsi(i)
       flu=0
       if(izeta(i).ne.1)then
        alil=-Phig*mass(i)*(1+ali)/temper(i)
        flu=densit(i)*sqrt( temper(i)/mass(i) )
        if(ali.gt.-1)flu=flu*(1+alil)*exp(-alil)
        if(iz(i).lt.0)flu=-flu
        sflux=sflux+flu
       endif
       print *,i,' alpha= ',ali,' flux= ',flu,' fnorm=',fnorm(i)
      enddo
      print *,'                   somme flux= ',sflux
c Phi neutre de 1000 en 1000 km    et h limite oxygene
      PhiNr=Phi1
      PhiN(1)=Phi1
      do ih=2,70
        h0=ih*1e6
c  Newton Raphson for neutral value at h0
        resid(1)=0
        do ip=1,9999
         resid(2)=resid(1)
         PhiNr=PhiNr-0.01
         resid(1)=0
         ti=0
         do j=1,izmax
            call density(h0,hb,den,nom(j),densit(j),fnorm(j),dden,
     &      PhiNr,Phi1,
     &      mass(j),temper(j),iz(j),ic(j),ipsi(j),ixi(j),izeta(j),0)
c          oxygene!!!
             if( abs(mass(j)/mp-16).lt.0.001 )then
              if(den.gt.denox)then
               den=denox
               dden=0
              endif
             endif
          resid(1)=resid(1)-ee*iz(j)*den
          ti=ti-ee*iz(j)*dden
         enddo
         if(resid(1).lt.0)then
          dres=   (resid(1)-resid(2))/(-0.01)
          PhiNr=PhiNr-resid(1)/dres
          resid(1)=0
          ti=0
          do j=1,izmax
            call density(h0,hb,den,nom(j),densit(j),fnorm(j),dden,
     &      PhiNr,Phi1,
     &      mass(j),temper(j),iz(j),ic(j),ipsi(j),ixi(j),izeta(j),0)
c           oxygene!!!
             if( abs(mass(j)/mp-16).lt.0.001 )then
              if(den.lt.denox)then
               denox=den
               hox=h0
              else
               den=denox
               dden=0
              endif
             endif
           resid(1)=resid(1)-ee*iz(j)*den
           ti=ti-ee*iz(j)*dden
          enddo
         goto 4
         endif
        enddo
   4   continue
c      print '('' neutral Phi at '',1p,e8.1,'' = '',0p,f9.5,''V , err=
c    &'',1p,E11.4,'' derr= '',e11.4,'' den O+ ='',E11.4 )',
c    &  h0,PhiNr,resid(1),dres,denox
       PhiN(ih)=PhiNr
       if(PhiN(ih).lt.PhiN(ih-1)-5)hcrit=h0-1e6
      enddo
c Phi neutre de 10 en 10 km
      ihx=hcrit/1e6-0.1
      PhiNr=PhiN(ihx)
      Phi0n=PhiNr
      do ih=1,101
        h0=hcrit+(ih-1)*1e4
        call pneutr(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        PhiNr,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
c      print '('' neutral Phi at '',f9.0,'' = '',f9.5,''V , err=
c    &'',1p,E11.4,'' derr= '',e11.4,'' den O+ ='',E11.4 )',
c    &  h0,PhiNr,err,derr,denox
       if(PhiNr.lt.Phi0N-5)then
           hcrit=h0-1e4
           PhiNn=Phi0N
       endif
      Phi0N=PhiNr
      enddo
c Phi neutre de 100 en 100 m
      ihx=hcrit/1e6-0.1
      PhiNr=PhiN(ihx)
      Phi0n=PhiNr
      do ih=1,101
        h0=hcrit+(ih-1)*100
        call pneutr(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        PhiNr,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
c      print '('' neutral Phi at '',f9.0,'' = '',f9.5,''V , err=
c    &'',1p,E11.4,'' derr= '',e11.4,'' den O+ ='',E11.4 )',
c    &  h0,PhiNr,aux1,dres,denox
       if(PhiNr.lt.Phi0N-5)then
           hcrit=h0-100
           PhiNn=Phi0N
       endif
      Phi0N=PhiNr
      enddo
c charge pres de hcrit: Taylor
      Phi0=PhiNn
      ipn=2*Phi0
         aux= charge(hcrit,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
      tayl(1)=aux
      print *,' charge en ',Phi0,'V ,   hcrit= ',hcrit,' : ',aux
         aux1= charge(hcrit+100,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
         aux2= charge(hcrit-100,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
      tayl(2)=(aux1-aux2)/200
      print *,' d/dh = ',tayl(2)
      tayl(3)=(aux1-2*aux+aux2)/20000
      print *,' 0.5 d2/dh2 = ',tayl(3)
         aux1= charge(hcrit,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0+0.1,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
         aux2= charge(hcrit,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0-0.1,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
      tayl(4)=(aux1-aux2)/0.2
      print *,' d/dPhi = ',tayl(4)
      tayl(5)=(aux1-2*aux+aux2)/0.02
      print *,' 0.5 d2/dPhi2 = ',tayl(5)
         aux1= charge(hcrit+100,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0+0.1,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
      tayl(6)=(tayl(1)-tayl(2)*100-tayl(3)*10000-tayl(4)*0.1
     &    -tayl(5)*0.0001)/10
      print *,'  d2/dhdPhi = ',tayl(6)
      aux=tayl(1)-tayl(4)*tayl(4)/(4*tayl(5))
      aux1=tayl(2)-tayl(4)*tayl(6)/(4*tayl(5))
      aux2=tayl(3)-tayl(6)*tayl(6)/(4*tayl(5))
      print *,' trinome:',aux,aux1,aux2
      aux2=hcrit+( -aux1+sign(one,aux)*sqrt( aux1*aux1-aux*aux1))/aux
      print *,' h critique= ',aux2
c     print  '(7x,8(f9.0))',( (hcrit+1000*(ih-4)),ih=1,8)
      print *,'  fichier postscript? 1=oui '
      read *,ipos
      if(ipos.ne.1)goto 42
      ipn=2*Phi1
      do ip=1,ipn-1
       Phi0=0.5*(ipn-1-ip)
       iPhi0=Phi0
       if(abs(Phi0-iPhi0).lt.0.01)
     &   write(12, '(''0 0 0 setrgbcolor '',f6.1,'' '',f6.1,''  m ('',
     &            f6.1,'' ) show '')')
     &    20.0,100+20*Phi0-3,Phi0
       do ih=1,100
c        h0=hcrit+(ih-1)*10000
         h0=0.5*ih*1E6
         if(h0.le.hb)goto 405
         aux= charge(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
c        xi(ih)=aux
        if(aux.gt.0)then
         write(12, '(''1 0 0 setrgbcolor '',f6.1,'' '',f6.1,''  m  '',
     &            f6.1,'' '',f6.1,''  l  '',
     &            f6.1,'' '',f6.1,''  l  '',
     &            f6.1,'' '',f6.1,''  l  '',
     &            f6.1,'' '',f6.1,''  l stroke '')')
     &    50+h0*1e-5-2,100+20*Phi0,50+h0*1e-5+2,100+20*Phi0,
     &    50+h0*1e-5,100+20*Phi0,50+h0*1e-5,100+20*Phi0+2,
     &    50+h0*1e-5,100+20*Phi0-2
        else
       write(12, '(''0 0 1 setrgbcolor  '',f6.1,'' '',f6.1,''  m  '',
     &            f6.1,'' '',f6.1,''  l stroke '')')
     &    50+h0*1e-5-2,100+20*Phi0,50+h0*1e-5+2,100+20*Phi0
        endif
  405  continue
       enddo
c      print '(2x,f5.2,1p,8(e9.2))',Phi0,( xi(ih),ih=1,8)
      enddo
      do ih=1,6
       h0=1e7*(ih-1)
       write(12, '(''0 0 0 setrgbcolor  '',f6.1,'' '',f6.1,''  m  '',
     &            f6.1,'' '',f6.1,''  l stroke '')')
     &    50+h0*1e-5,95.0,50+h0*1e-5,90.0
      enddo
      write(12,'(''showpage'')')
   42 print *,' calcul charge,  read h, Phi   stop si h<0'
       read *,h0,Phi0
       if(h0.lt.0)goto 44
         aux= charge(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
       print *,' h= ',h0,' Phi= ',Phi0,' charge= ',aux
      goto 42
   44 print *,' calcul de Phi neutre en d''autres h: read h stop si<0'
       read *,h0
       if(h0.lt.0)goto 5
       i0=h0/1E6-0.1
       Phi0=PhiN(i0)+0.05
        call pneutr(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
       print '('' neutral Phi at '',1p,e12.5,'' = '',0p,f9.5,''V , err=
     &'',1p,E11.4,'' derr= '',e11.4)',h0,Phi0,err,derr
  441 print *,' test balistique: read Phi'' , retour neutre si <-999'
      read *,aux1
      if(aux1.lt.-999)goto 44
      PhiNn=Phi0
       do i=1,100
         aux2= charge(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        PhiNn,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
c derivee seconde a partir du laplacien
         aux2=aux2-aux1/(h0+RE)
         PhiNn=PhiNn+aux1+aux2/2
         aux1=aux1+aux2
c         if(mod(i,10).eq.0)print *,h0,PhiNn,aux1,aux2
         if(mod(i,10).eq.0)
     &   print '(a1,''psdots[dotstyle=.,linecolor=black,dotsize=2pt]('',
     &    f6.3,'','',f6.3,'')'')',char(92),(h0-2.3E7)*1E-5,0.8*PhiNn
         h0=h0+1
       enddo
      goto 441
    5 print *,' Poisson: h0,hN,Npo? (stop si <=0)'
      read *,h0,hN,Npo
      if(Npo.le.0)stop
       i0=h0/1E6-0.1
       Phi0=PhiN(i0)+0.05
        call pneutr(h0,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi0,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
       iN=hN/1E6-0.1
       PhiNn=PhiN(i0)+0.05
        call pneutr(hN,hb,hox,denox,nom,densit,fnorm,izmax,
     &        PhiNn,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
      print *,' h0=',h0,' hN=',hN,' Phi0=',Phi0,'PhiN=',PhiNn
      do i=0,Npo
       xi(i)=(h0+hN-(hN-h0)*cos(i*pi/Npo)  )/2
      enddo
         if(Npo.lt.8)print '(1p,7e10.3)',(xi(j),j=0,Npo)
         print *,' '
       test=0
      do i=1,Npo-1
        call weights1(xi(i),xi,Npo,2,fornb,ndim)
         if(Npo.lt.8)print '(1p,7e10.3)',(fornb(j,2),j=0,Npo)
        ti=0
        do j=1,Npo-1
         D2(i,j)=fornb(j,2)+(2/(RE+xi(i)))*fornb(j,1)
         ti=ti+D2(i,j)*(xi(j)-h0)*(xi(j)-hN)
        enddo
       test=test+abs( ti-2 -2*(2*xi(i)-h0-hN)/(xi(i)+RE) )
      enddo
      print *,' test Laplace: ',test
c  N.B.   D2 donne la derivee 2de aux points interieurs d'une
c           fonction nulle aux extremites!
c  first approximation for Phi
      do i=0,Npo
c      i0=xi(i)/1E6-0.1
c      Phi(i)=PhiN(i0)+0.05
c       call pneutr(xi(i),hb,hox,denox,nom,densit,fnorm,izmax,
c    &        Phi(i),Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
       Phi(i)=Phi0+(xi(i)-h0)*(PhiNn-Phi0)/(hN-h0)
      enddo
         if(Npo.lt.8)print '(1p,7e10.3)',(Phi(j),j=0,Npo)
      itertot=0
1     print *,' how many further iterations? 0=stop, relax<=1'
      read *,iter,relax
      if(iter.le.0)goto 2
      do iterm=1,iter
       itertot=itertot+1
c residue F-Phi''  and A=dF/dPhi
       Fmax=0
       rmax=0
       do i=1,Npo-1
         aux= charge(xi(i),hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi(i)-0.01,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
         aux1= charge(xi(i),hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi(i),Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
        Fmax=max(Fmax,abs(aux1))
        do j=1,Npo-1
c  radial part of Laplace operator:
         mat(i,j)=D2(i,j)
         resid(i)=aux1-D2(i,j)*
     &    ( Phi(j)-  Phi0-(PhiNn-Phi0)*(RE+hN)*(xi(j)-h0)/
     &          ((hN-h0)*(xi(j)+RE))  )
        enddo
       rmax=max(rmax,abs(resid(i)))
c         do j=1,izmax
          mat(i,i)=mat(i,i)-(aux1-aux)/0.01
c         enddo
       enddo
       print *,itertot,' Laplace: ',Fmax,' err: ',rmax
       call GAUSS(mat,Phicor,resid,Npo-1,ndim)
       do i=1,Npo-1
        Phi(i)=Phi(i)+relax*Phicor(i-1)
       enddo
         if(Npo.lt.8)print '(1p,7e10.3)',(Phi(j),j=0,Npo)
      enddo
      goto 1
c  some values of Phi and Phi'
2     print *,' points du graphe? entre 10 et 500'
      read *,ngra
      if(ngra.lt.10.or.ngra.gt.500)stop
      scalx(1)=-990
      scalx(2)=-990
      scalx(3)=-990
      scalm(1)=990
      scalm(2)=990
      scalm(3)=990
      fgra=1.0/ngra
      do i=5,ngra-5
       plot(i,1)=0
       plot(i,2)=0
       plot(i,3)=0
       hi=h0+fgra*i*(hN-h0)
       call weights1(hi,xi,Npo,2,fornb,ndim)
       do j=0,Npo
        plot(i,1)=plot(i,1)+fornb(j,0)*Phi(j)
        plot(i,2)=plot(i,2)-fornb(j,1)*Phi(j)
        plot(i,3)=plot(i,3)+(fornb(j,2)+2*fornb(j,1)/(hi+RE))*
     &    ( Phi(j)-  Phi0-(PhiNn-Phi0)*(RE+hN)*(xi(j)-h0)/
     &          ((hN-h0)*(xi(j)+RE))  )
       enddo
       scalx(1)=max(scalx(1),plot(i,1))
       scalx(2)=max(scalx(2),plot(i,2))
       scalx(3)=max(scalx(3),plot(i,3))
       scalm(1)=min(scalm(1),plot(i,1))
       scalm(2)=min(scalm(2),plot(i,2))
       scalm(3)=min(scalm(3),plot(i,3))
      enddo
      print *,' min Phi = ',scalm(1), ' max= ',scalx(1)
      print *,' min -dPhi/dh = ',scalm(2), ' max= ',scalx(2)
      print *,' min Laplacien Phi = ',scalm(3),' max= ',scalx(3)
      print *,char(92),'psset{unit=1cm}'
      do j=1,3
       scal=scalx(j)-scalm(j)
       print *,char(92),'begin{pspicture}(0,-0.5)(15,5)'
       print '(a1,''psline'',4(''('',f6.2,'','',f6.2,'')''))',char(92),
     &0,4*(plot(5,j)-scalm(j))/scal,0,-4*scalm(j)/scal,15.0,-4*scalm(j)
     &/scal,15.0,4*(plot(ngra-5,j)-scalm(j))/scal
       print '(a1,''uput[45](0,4){$'',1pe11.4,''$}'')',char(92),scalx(j)
       print '(a1,''uput[0](0,0){$'',1pe11.4,''$}'')',char(92),scalm(j)
       print *,char(92),'pscurve[linewidth=0.01]%'
       km=4
       do i=0,ngra/5
        if(i.eq.ngra/5)km=0
        print '(5(''('',f6.2,'','',f6.2,'')''),''%'')',(0.03*(5*i+k),
     &   4*(plot(5*i+k,j)-scalm(j))/scal,k=0,km)
       enddo
       print *,char(92),'pscurve[linewidth=0.05]%'
       do i=0,Npo
        hi=15*(xi(i)-h0)/(hN-h0)
        k=hi/0.03
        print '(''('',f6.2,'','',f6.2,'')%'')',hi,
     &   4*(plot(k,j)-scalm(j))/scal
       enddo
       do i=0,Npo
        hi=15*(xi(i)-h0)/(hN-h0)
        k=hi/0.03
        print '(a1,''psdots[dotsize=0.12,dotstyle=o]('',f6.2,'','',
     &    f6.2,'')'')',char(92),hi,
     &   4*(plot(k,j)-scalm(j))/scal
       enddo
       print *,char(92),'end{pspicture}'
      enddo
      stop
      end
c-----------------------------------------------------------------
      BLOCK DATA
c constantes physiques
        real e,ee,Tk,mp,me,RE
        COMMON /cphys/e,ee,Tk,mp,me,RE
c e=el charge, ee=e/8.85E-12, Tk=Boltz., mp=mass p, me=mass e=mp/1836
        data e,ee,Tk,mp,me,RE/1.602E-19,1.810E-08,1.38E-23,1.67E-27,
     &   9.096E-31,6.371E6/
      END
c-----------------------------------------------------------------
      subroutine density(h,hb,den,nom,denb,fn,dden,Phi,Phi1,mass,temp,
     &        iz,ic,ipsi,ix,izet,iverb)
c  mass, temp, iz,ic,ipsi,ix,izet= mass etc.  of particle
c   (temp = k temperature )
c  hb,denb = given reference altitude and density
c  h = given altitude
c den=sought density     dden= derivative wrt Phi
      real*8 h,hb,den,denb,fn,dden,Phi,Phi1,Phig,mass,temp,
     &     y,q,one,pi,aux,aux1,eta,fk2,vinf2,xsq
      character*4 nom
        real e,ee,Tk,mp,me,RE
      COMMON /cphys/e,ee,Tk,mp,me,RE
      one=1
      pi=4*atan(one)
      y=(RE+hb)/(RE+h)
      Phig=-9.81*RE/(1+hb/RE)
      q=(mass*Phig*(y-1)+iz*e*(Phi-Phi1))/temp
c  barometric case
        den=denb* exp(-q)
        dden=-e*iz*den/temp
      if( ((ic.ne.1) .or. (abs(y-one).lt.0.00000001)).and.iverb.gt.0 )
     &then
        print '(a4,'' q='',1p,e12.5,'' ic='',i2,'' den='',e12.5,
     &'' dden='',e12.5)',nom,q,ic,den,dden
        return
      endif
c  exospheric case
      eta=y**3
      if(q.ge.0)then
       aux=(-mass*Phig*y-iz*e*Phi)/temp
       vinf2=aux
c      if(aux.lt.0)print *,' ! Vinf2=',aux
       if(aux.lt.0)aux=-mass*Phig/(temp*1e10)
       aux1=( aux + eta*(mass*Phig+iz*e*Phi1)/temp )/(one-eta)
       xsq=aux1
c      if(aux1.lt.0)print *,' ! X2=',aux1
       if(aux1.lt.0)aux1=-mass*Phig/(temp*1e10)
       den=den*( ipsi+(izet-ipsi)*fk2(sqrt(aux))
     &   -sqrt(one-eta)*( ipsi-ix-(ipsi+izet-2*ix)*fk2(sqrt(aux1)))
     &         *exp(-eta*q/(one-eta))    )/fn
        if(iverb.gt.0)
     &  print '(a4,'' q='',1p,e12.5,'' Vinf2='',e12.5,'' xsq2='',e12.5,
     &'' den='',e12.5)',nom,q,vinf2,xsq,den
      else
       den=den*( izet+(ipsi-izet)*(   fk2(sqrt(-q))
     &              -sqrt(one-eta)*exp(-eta*q/(one-eta))
     &                              *fk2(sqrt(-q/(one-eta)))
     &                             )   )/fn
      if(iverb.gt.0)
     & print '(a4,'' q='',1p,e12.5,''       '',12x,''      '',12x,
     &'' den='',e12.5)',nom,q,den
      endif
      return
      end
c-----------------------------------------------------------------
      real function charge*8(h,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
      parameter (ndim=75,ndat=12)
      real*8 Phi,hb,Phi1,densit(ndat),mass(ndat),temper(ndat),aux,aux1,
     &  h,fnorm(ndat),denox,hox,den,dden
      integer iz(ndat),ic(ndat),ipsi(ndat),ixi(ndat),izeta(ndat)
      character*4 nom(ndat)
        real e,ee,Tk,mp,me,RE
      COMMON /cphys/e,ee,Tk,mp,me,RE
        aux=0
         do j=1,izmax
            call density(h,hb,den,nom(j),densit(j),fnorm(j),dden,
     &      Phi,Phi1,
     &      mass(j),temper(j),iz(j),ic(j),ipsi(j),ixi(j),izeta(j),0)
c          oxygene!!!
             if( (abs(mass(j)/mp-16).lt.0.01).and.(h.gt.hox) )then
               den=denox
               dden=0
             endif
          aux=aux-ee*iz(j)*den
         enddo
        charge=aux
      return
      end
c-----------------------------------------------------------------
      subroutine pneutr(h,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta,err,derr)
c remplace Phi par Phi neutre (en decroissant)
      parameter (ndim=75,ndat=12)
      real*8 Phi,hb,Phi1,densit(ndat),mass(ndat),temper(ndat),aux,aux1,
     &  h,fnorm(ndat),denox,hox,den,dden,err,derr,charge
      integer iz(ndat),ic(ndat),ipsi(ndat),ixi(ndat),izeta(ndat)
      character*4 nom(ndat)
        real e,ee,Tk,mp,me,RE
      COMMON /cphys/e,ee,Tk,mp,me,RE
         err= charge(h,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
        if(err.lt.0)then
         print *,' Pneutr: Phi=',Phi,' trop petit ',err
         return
        endif
        do ip=1,9999
         aux1=err
         Phi=Phi-0.01
         err= charge(h,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
         if(err.lt.0)then
          derr=   (err-aux1)/(-0.01)
          Phi=Phi-err/derr
         err= charge(h,hb,hox,denox,nom,densit,fnorm,izmax,
     &        Phi,Phi1,mass,temper,iz,ic,ipsi,ixi,izeta)
         goto 45
         endif
        enddo
  45   continue
       return
       end
c-----------------------------------------------------------------
      real function fk2*8(x)
c  1/2-K2 !     =  (  x exp(-x2) +int de x a infty de exp(-t2)dt )
c                            /sqrt(pi)
c                       large x:     sim (x+1/(2x))exp(-x2)/spi
       real*8 x,one,spi
       one=1
       spi=sqrt( 4*atan(one) )
       fk2=erfc(x)/2 +x*exp(-x*x)/spi
      return
      end
c-----------------------------------------------------------------
      subroutine weights1(xi,x,n,m,c,ndim)
c
c  Cf. B. Fornberg, A Practical Guide to Pseudospectral Methods,
c   Cambridge University Press, 1996,  paperback: 1999.
c
c  B. Fornberg, Calculation of weights in finite difference
c     formulas, SIAM Rev. vol. 40 no 3 , 685-691, Sept. 1998.
c
c+-------------------------------------------------------------------+
c|  INPUT PARAMETERS:                                                |
c|    XI  POINT AT WHICH THE APPROXIMATIONS ARE TO BE ACCURATE       |
c|    X   X-COORDINATES FOR GRID POINTS, ARRAY DIMENSIONED X(0:N)    |
c|    N   THE GRID POINTS ARE AT X(0),X(1),...,X(N) (I.E. N+1 IN ALL)|
c|    M   HIGHEST ORDER OF DERIVATIVE TO BE APPROXIMATED             |
c|    ndim    see dimensions of c                                    |
c|                                                                   |
c| OUTPUT PARAMETER:                                                 |
c|    C   WEIGHTS, ARRAY DIMENSIONED C(0:Ndim,0:M).                     |
c|        ON RETURN, THE ELEMENT C(J,K) CONTAINS THE WEIGHT TO BE    |
c|        APPLIED AT X(J) WHEN THE K:TH DERIVATIVE IS APPROXIMATED   |
c|        BY A STENCIL EXTENDING OVER X(0),X(1),...,X(N).            |
c+-------------------------------------------------------------------+
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(0:N),C(0:Ndim,0:M)
      C1=1
      C4=X(0)-XI
      DO 10 K=0,M
         DO 10 J=0,N
   10       C(J,K)=0
      C(0,0)=1
      DO 50 I=1,N
         MN=MIN(I,M)
         C2=1
         C5=C4
         C4=X(I)-XI
         DO 40 J=0,I-1
            C3=X(I)-X(J)
            C2=C2*C3
         DO 20 K=MN,1,-1
   20       C(I,K)=C1*(K*C(I-1,K-1)-C5*C(I-1,K))/C2
         C(I,0)   =-C1*C5*C(I-1,0)/C2
            DO 30 K=MN,1,-1
   30          C(J,K)= (C4*C(J,K)-K*C(J,K-1))/C3
   40       C(J,0)=C4*C(J,0)/C3
   50    C1=C2
      RETURN
      END
c-----------------------------------------------------------------
      SUBROUTINE GAUSS(AM,XV,BV,ND,ndim)
C RESOLUTION DE AX=B PIVOTAGE PARTIEL:PAX=PB
        REAL*8 AM(ndim,50),XV(50),BV(50),B1,B2,PMAX,CCABS
      DIMENSION IP(200)
        ISIGN=1
      N1=ND-1
      DO 10 I=1,ND
10    IP(I)=I
      DO 1 J=1,ND
      PMAX=0
      DO 11 I=J,ND
       CCABS=ABS(AM(IP(I),J))
      IF(CCABS.LE.PMAX)GOTO 11
      PMAX=CCABS
      IPP=I
11    CONTINUE
      IPJ=IP(J)
        KP=IP(IPP)-IPJ
        IF(KP.NE.0)ISIGN=-ISIGN
      IP(J)=IP(IPP)
      IP(IPP)=IPJ
      B1=1/AM(IP(J),J)
      JJ=J+1
      IF(N1)12,12,13
12    XV(1)=B1*BV(1)
      GOTO 99
13      IF(J.EQ.ND)GOTO 1
      DO 101 I=JJ,ND
      B2=AM(IP(I),J)*B1
      BV(IP(I))=BV(IP(I))-B2*BV(IP(J))
      DO 101 K=JJ,ND
  101 AM(IP(I),K)=AM(IP(I),K)-B2*AM(IP(J),K)
1        CONTINUE
      XV(ND)=BV(IP(ND))/AM(IP(ND),ND)
      DO 2 I=1,N1
      B1=BV(IP(ND-I))
      DO 3 J=1,I
3     B1=B1-AM(IP(ND-I),ND-I+J)*XV(ND-I+J)
2     XV(ND-I)=B1/AM(IP(ND-I),ND-I)
99    RETURN
      END
