Quantenmechanik/ Boson-Fermion

Identische Teilchen Bearbeiten

Ein System aus mehreren Teilchen hat quantenmechanisch nicht etwa mehrere Wellenfunktionen, sondern eine gemeinsame im Tensorprodukt der Einteilchen-Hilbert-Räume. Nur wenn für Teilchen A bei x und B bei y keine Wechselwirkung   im Hamilton-Operator auftritt, ist ein Separations-Ansatz   erfolgreich. Allgemein gibt es Linearkombinationen von solchen Produkten, also verschränkte Teilchen:  

Teilchen wie Protonen, Elektronen oder ganze Atome haben nicht das geringste Merkmal, mit dem man sie individuell verfolgen kann. Gleiche Teilchen sind ununterscheidbar. Sie haben keine Pfade wie in der klassischen Physik. Sind Teilchen A und B identisch wie die zwei Elektronen des Helium-Atoms, dann stellt sich schon ohne Wechselwirkung eine Frage.

Die Wellen   und   sind degeneriert, es könnte also sein, dass ein ganzes Kontinuum von Funktionen

 

die selbe Physik beschreibt. Der Permutations-Operator   hätte dann keine definierten Eigenwerte?

Oder aber P hat eine einfache Darstellung,  , das heißt die Zweiteilchen-Funktion ist entweder symmetrisch oder antisymmetrisch beim Vertausch der Koordinaten, also mit Eigenwert 1 oder -1 zu P. Auf jeden Fall hätte P einen Eigenwert auf dem Einheitskreis, so dass die observablen Betragsquadrate invariant sind beim Paarvertausch?

Die Antwort der Natur ist kategorisch und wurde von Pauli formuliert. Die Mehrteilchen-Zustände N identischer Teilchen mit halbzahligem Spin, wie der Elektronen, haben als Darstellung der Permutationsgruppe den Faktor -1 bei Paarvertauschung. Teilchen mit dieser Antisymmetrie sind Fermionen. Teilchen mit ganzzahligem Spin, die Bosonen, haben den P-Eigenwert 1, sie haben also vertauschungs-invariante N-Teilchen-Zustände.

Unmittelbare Konsequenz ist das Pauli-Verbot: keine zwei Fermionen passen auf den selben Quantenzustand:  . Mittelbare Konsequenz ist der Aufbau der Atome der Kernladung Z in Schalen. Die Elektronen füllen nacheinander die aufsteigende Reihe der Energie-Niveaus, die bei der Berechnung eines Elektrons im Coulomb-Potenzial herauskommen.

Die Ursache des Spin-Statistik-Prinzips ist in der relativistischen Quantenfeldtheorie zu suchen. Die identischen Teilchen sind alle die elementaren Anregungen des selben Quantenfeldes. Der Spin ergibt sich aus einer Klassifizierung der Darstellungen der Lorentz-Gruppe. Die Vielteilchen-Zustände sollen einen Hilbert-Raum aufspannen (positives Skalarprodukt) und einen stabilen Grundzustand haben (beschränktes Spektrum der Energie). Nur die experimentell untermauerte Korrespondenz von Spin und Vertauschungs-Symmetrie ist dann widerspruchsfrei möglich.

Statistik Bearbeiten

Angenommen, es gebe   Plätze oder Einzelniveaus für einzelne Teilchen. Jedem Teilchen gehöre ein Vektorraum der Dimension I mit einer Basis von Einteilchen-Zuständen  .

Naiv besteht ein Basissystem der Zustände von   verfolgbaren Teilchen aus allen Abbildungen  . Es gibt   Vektoren in der Basis des gewöhnlichen Tensorproduktes  . Sie beschreibt das System:

  erzeugt  .

Die Dimension ist die Zahl der linear unabhängigen Zuständsvektoren.

Mit Fermionen sind für jedes i nur zwei Besetzungen erlaubt, 0 oder 1. Falls  , sind die unabhängigen Zustände die Menge aller Teilmengen der Größe N, insgesamt   Zustände. Soviel Möglichkeiten gibt es, N Plätze aus den I verfügbaren auszuwählen. Der Vektorraum des Systems wird aufgespannt von allen total antisymmetrischen Tensorprodukten, ist also ein Teilraum  , definiert wie folgt.

Sei   die Gruppe aller N! Permutationen von  . Die Gruppe hat zwei elementare Homomorphismen zu den ganzen Zahlen,

  • den trivialen,   und
  • das Permutations-Vorzeichen   genau dann wenn
p aus einer ungeraden Zahl von Nachbar-Vertauschungen gebaut ist, 1 sonst.

Alle Paarvertauschungen, nicht nur die Transpositionen von Nachbarn, haben negatives  . Eine beliebige Permutation kann in ein Produkt aus Zyklen zerlegt werden. Ein Zyklus der Größe   hat das Vorzeichen  .

Der Antisymmetrie-Operator A auf   ist die lineare Fortsetzung von:

 

Der Vektorraum des Fermionsystems ist die Bildmenge  . A ist eine Projektion,  . Als Projektion hat A die Tendenz, Vektoren zu "verkürzen" und erhält daher nicht die Normen.

In der Ortsdarstellung geht die Vielteilchen-Welle aus Produkten hervor. Man wende den Operator A an auf  .

 

Eine für orthogonale Produkte norm-erhaltende Variante hat folgende Form:

 .

Diese Schreibweise einer N-Fermion-Wellenfunktion ist bekannt under dem Namen Slater-Determinante. Alle beteiligten Einteilchen-Wellen müssen linear unabhängig sein, damit sie nicht verschwindet. Die Slater-Determinanten eines Einteilchen-Orthonormalsystems bilden ein Vielteilchen-Orthonormalsystem.

Die Argumente   können allgemeiner Tupel sein aus Ortskoordinaten, Spin-Indizes und anderen Quantenzahlen. Bei N Elektronen können eventuell nur die Spin-Faktoren der Wellenfunktion antisymmetrisch sein; dann werden die Orts-Faktoren symmetrisch für die geforderte totale Antisymmetrie. Alle möglichen Kombinationen der Symmetrien mit Tupeln können mit der Darstellungstheorie der Permutationsgruppen (Young-Tableaus etc) klassifiziert werden.

Bei Bosonen darf man N "Ziehungen mit Zurücklegen" aus den I Einzelniveaus machen, was der Zahl von Möglichkeiten   entspricht.

Die immer wieder erquickliche Umformulierung beweist diese Kombinatorik:

Die N Kugeln sollen auf I Fächer   verteilt werden. Die Fächer kriegen (I-1) Trennwände, so dass jede Verteilung genau einer Folge aus (N+I-1) Dingen entspricht, entweder Kugel oder Wand. Gesucht ist die Zahl aller Möglichkeiten, die (I-1) Wände zu setzen. Das ergibt  .

Der Boson-Vektorraum ist der Unterraum der total symmetrischen Tensorprodukte in  , er ist die Bildmenge des Symmetrierung-Operators S. Operator   ist definiert wie folgt (linear zu erweitern):

 

Auch der Operator S ist eine Projektion. Eine norm-erhaltende Version für Produkte orthogonaler Funktionen wäre:  .

Der Grundzustand mit vielen Bosonen ist der, wo alle auf dem selben niedrigsten Energie-Niveau hocken: Bose-Einstein-Kondensation. Er hat mehr Gewicht bei insgesamt   im Vergleich zu den   klassisch markierbaren möglichen Zustands-Dimensionen.

Die Multi-Fermion- und Multi-Boson-Zustände liegen in Eigenräumen zum Eigenwert 1 der Operatoren A bzw. S. Der Hamilton-Operator und alle Operatoren von Observablen sollten diesen Zustandsraum nicht verlassen können und sollten mit allen Permutations-Operatoren von Teilchenpaaren (Transpositionen) kommutieren. Es ist auch nicht möglich, dass ein System aus gleichen Teilchen sein Verhalten dynamisch ändert, von Symmetrie zu Antisymmetrie oder zurück.

(In der spekulativen Physik gibt es Super-Symmetrien mit Operatoren, die Bosonen und Fermionen ineinander umwandeln! Bisher wurde nichts dergleichen in der Natur gefunden.)

Auf der direkten Summe aller Boson-Vektorräume zu allen Teilchenzahlen

 

gibt es ein Produkt. Es wird auf Monomen von Produktzuständen von N bzw. M Teilchen als   definiert, bilinear fortzusetzen auf alle Polynome. Das Produkt ist assoziativ und kommutativ; der Boson-Raum ist damit eine Kommutative Algebra. Die erste Dimension von komplexen Zahlen wird hinzugenommen und ist physikalisch der Zustand des Vakuums, das null Teilchen enthält. Damit hat die Algebra ein Einselement.

Auf der direkten Summe der Fermion-Vektorräume zu allen Teilchenzahlen

 

ist das entsprechende Produkt auf Monomen von N bzw. M Teilchen als   definiert, bilinear fortzusetzen auf alle Polynome. Das Produkt ist assoziativ und alternierend. Der Fermion-Raum ist damit eine Grassmann-Algebra. Die Vertauschung funktioniert mit der Vorzeichenregel:

 .

Zwei Monome mit ungeraden Einteilchen-Faktoren antikommutieren, alle anderen Paare kommutieren.

Die beiden Algebren haben einen Grad bei Monomen, die Zahl der elementaren Faktoren (wie bei gewöhnlichen Polynomen der Homogenitäts-Grad). Sie werden graduierte Algebren genannt. Der Grad ist additiv bei Multiplikation von homogenen Elementen. Viele Operatoren auf den Algebren variieren den Grad um 0,+1,-1 Einheit. Manche agieren auf Produkten mit einer Leibniz-Regel, ähnlich zur Ableitung vom Produkt zweier Funktionen; sie heißen dann Derivationen bzw. Anti-Derivationen.

Exponential-Verteilung Bearbeiten

In den einfachsten klassischen Modellen der statistischen Mechanik sind die Energie-Zustände der Teilchen exponentiell verteilt. Die relative Häufigkeit der Teilchen mit Energie  .

Quantenmechanisch gibt es weit weniger 'Freiheitsgrade', wie schon aus der Abzählung der Basen hervorging. Pauli verbietet den Fermionen die Doppelbelegung eines Niveaus, und auch bei Bosonen zählen gleiche Partikelbelegungen der Niveaus nur als ein einziger Zustand.

Ein Multi-Boson-Zustand ist eine Folge   ganzer Zahlen; seine Häufigkeit ist  . Interessant ist die relative Häufigkeit von Niveau  , wofür der Erwartungswert von   aus der Verteilung   zu berechnen ist.

 

Faktoren   können aus allen   gestrichen werden. Es bleibt mit  :

 .
 
 

Ein Multi-Fermion-Zustand ist eine Folge   die nur aus den Werten 0 und 1 besteht. Die Wahrscheinlichkeit für ein Teilchen auf dem Niveau i wird dann ganz analog mit Summen über {0,1} statt   :

 

??? Normierung? Diese Verteilung wird nahezu Null für E_i weit oberhalb der thermischen Energie k_BT=1/\beta, nahezu konstant für Energien weit darunter? negative?

Zusammenfassung

  • Klassische Verteilung  .
  • Fermi-Dirac-Verteilung  
  • Bose-Einstein-Verteilung  

Periodensystem Bearbeiten

 
Radiale Dichte der Elektronenladung von Edelgas-Atomen. Zeigt den Aufbau des Elektronensystems in Schalen

Nach dem Pauli-Prinzip bauen sich die Atome um einen Kern der Ladung Z auf. In erster Näherung belegt das erste Elektron das niedrigste Energieniveau von den Lösungen im Coulomb-Potenzial. Die weiteren Elektronen belegen höhere Plätze. Sie sehen ein verformtes Potenzial wegen der Abschirmung des Kerns durch alle anderen Elektronenwolken.

Die einfachste Näherung wäre, für alle Z Elektronen das selbe mittlere radiale Potenzial   einzusetzen. Der vereinfachte Hamilton-Operator wäre

 .

Die Wellenfunktionen kämen aus einer Basis der Slater-Determinanten von wasserstoff-ähnlichen Lösungen.

Mit besseren Methoden und Computer-Hilfe konnte die Wellenmechanik nach Schrödinger und Pauli die Grundzustände aller Atome des Periodensystems naturgetreu berechnen. Die Entartung aller (n,l) zu gleichem n wird durch das kollektive Potenzial aufgehoben. Die Reihenfolge, in der die Paare (n,l) das höchste belegte Niveau stellen, konnte genau bestätigt werden.

Die folgende Tabelle zeigt die Ordnung der Niveaus mit Hauptquantenzahl 1-6 und Drehimpuls 0-3 (spdf), wie sie mit steigender Energie gefüllt werden. Bei gegebener Hauptquantenzahl haben die Drehimpulse s,p,d,f, zusammen mit den Paaren vom Elektronenspin, je 2,6,10,14 Zustände mit fast gleicher Energie. Jede Reihe beginnt mit einem Alkali-Atom und endet mit einem Edelgas. Vor dem Edelgas steht jeweils ein Halogen. Das Edelgas mit seiner stabilen Elektronenschale ist chemisch inert. Das Alkalimetall gibt reaktionsfreudig sein vereinzeltes Elektron ab. Genauso begierig nimmt das Halogen ein Elektron auf, um die Edelgaskonfiguration zu erreichen. Die theoretische und numerische Chemie ist nichts anderes als angewandte Quantenmechanik.

Niveaus Anzahl Elemente
1s 2 H He
2s 2p 2+6 Li ... F Ne
3s 3p 2+6 Na ... Cl Ar
4s 3d 4p 2+10+6 K ... Br Kr
5s 4d 5p 2+10+6 Rb ... I Xe
6s 4f 5d 6p 2+14+10+6 Cs ... At Rn

Die Konfigurationen der Elektronen der Alkalimetalle Natrium und Kalium und des Halogen-Atoms Chlor sehen schematisch aus wie folgt.

 
Konfiguration Na
 
Konfiguration Cl
 
Konfiguration K

Thomas-Fermi-Modell Bearbeiten

Dies ist ein vereinfachtes Modell für das kollektive Zentralpotenzial V(r), das von den Elektronen um einen Kern der Ladung Z gefühlt wird. Gefordert:

  • Asymptotisch für  
  • Für  
  •   (Elektron-Dichte im Grundzustand)
  •  

V ist am Zentrum von der Kernladung bestimmt und anderswo von der Ladungsdichte der gesamten Z Elektron-Wellen-Faktoren. Diese stecken in einer Slater-Determinante, die den Grundzustand annähert.

Nun kommt eine halb-klassische Idee hinzu: die Dichte von gebundenen Quantenzuständen im Phasenraum. In der WKB-Näherung sind die Energie-Niveaus des eindimensionalen Modells gequantelt nach der Regel:

 .

Die Ebene der Punkte   heißt der Phasenraum, seine allgemeine Version für N Freiheitsgrade hat die Dimension 2N.

Mit   kann die Quantelung so beschrieben werden: Das klassische System im Potenzialtopf oszilliert zwischen zwei Umkehrpunkten   und beschreibt geschlossene Kurven im Phasenraum  . Zuerst mit   für x von a bis b, dann zurück mit   für x von b nach a. Das Integral über eine Runde ist die umrundete Fläche in der  -Ebene.

 .

Einer Stufe der Energie E entspricht ein Zuwachs von "Volumen im Phasenraum"   vom Wert h. In drei Dimensionen wächst der Phasenraum um den Wert   zwischen Energieniveaus. Die Quantenmechanik schreibt also den gebundenen Fermionen wegen des Pauli-Verbots vor, mit einer konstanten Packungsdichte   den Phasenraum aufzufüllen, angefangen mit dem niedrigsten Energie-Niveau. Der Faktor 2 zählt die 2 Dimensionen des Spins.

Nach dem Thomas-Fermi-Modell sind die Energiezustände im 6-dimensionalen Phasenraum   lokalisiert mit der Dichte   falls die Energie   kleiner ist als das maximal besetzte Niveau  . Man darf   festlegen, denn eine Verschiebung ändert nicht das Verhalten von   am Nullpunkt. Der Faktor 2 berücksichtigt die zwei Spin-Richtungen der Elektronen. Die Elektronendichte im Raum:

 
 .

Die Dichte   hängt also von   ab und ist positiv wenn  , Null sonst. Es gilt also diese Differenzialgleichung für  :

 .

Es werden   durch dimensionslose Variablen   ersetzt:

 
  mit dem Grenzwert   wegen  .

Damit:   im Gebiet  .
Die Differenzialgleichung für u lautet in der Kugel wo   positiv ist:

 .

Wann kommt die erste Nullstelle   von u(x)? Ausgerechnet wird dazu eine Normierung von u, die daraus folgt dass   sich zu Z summiert:

 .

Daraus folgt   an der Nullstelle. Das weist darauf hin, dass   bei Unendlich liegt.   ist also zu lösen mit den zwei Randbedingungen  . Ein Skript im Anhang integriert die Gleichung numerisch und macht den Graphen als SVG-Datei.

Aus der Funktion u() folgen die radialen Potenziale   und Ladungsdichten  , die unter der Annahme   großer Atome angenähert beschreiben, wie dicht die Elektronen sich drängeln dürfen. Man definiert einen Atomradius   als die Kugel, in der der Bruchteil   der Elektronenladung ist. Typischerweise  , also die Kugel für alle bis auf ein Elektron. Gleichung:  . Mit dimensionslosem  , ergibt sich X als numerische Lösung von  . Ergebnis: R variiert recht langsam mit steigendem Z. Der Phasenraum wird schneller in der Impuls- als in der Orts-Dimension aufgeblasen.

Das Skript integriert numerisch die Gleichung   mit der Probeschuss-Technik wie folgt, ohne Anspruch auf hohe Präzision:

  • Intervall   wird in 100 gleiche Schritte d=0.01 geteilt.

Der typische Rechenschritt bestimmt   aus  

  • Intervall   wird als geometrische Folge in 1% große Schritte geteilt.

Hier ist u zur Funktion   umgeschrieben und die Rechenschritte sind die Differenz-Version von  .
Es wird   festgesetzt und   mit Dichotomie variiert, bis am Ende   herauskommt, von Rundungsfehlern abgesehen. Die SVG-Grafik zeigt die Funktion mit logarithmischer x-Achse.

Anhang: Skript Bearbeiten

Skript fürs Thomas-Fermi-Modell. Schreibt eine Datei thomferm.svg.

'''
Thomas-Fermi potential model.
solve x(f")^2 = f^3 on interval [0,inf), f(0)=1, f(1)=0
implemented direct shoot algorithm with log section:
f(x)=F(y); y=log(x); x=exp(y); f'=F'/x; f"= F"/x^2-F'/x^2=(F"-F')/(xx)
F"-F'= xx sqrt(F^3/x)= (F*x)^(3/2) = (F*exp(y))^(3/2)= z
(F(2)+F(0)-2F(1))/(dd)-(F(2)-F(0)/(2d) = z
(1-d/2)F(2)/dd = z + (2F(1)-F(0))/dd - F(0)/(2d) 
(1-d/2)F(2)= ddz+ 2F(1)- (1+d/2) F(0)
'''

from math import sqrt,exp,log
def div(p,q) : return int(p/q)

###### plot

class svgdump :
  def plotbegin(sd, w,h,f) : # args: width height zoomfactor
    W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
    s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
     '<!DOCTYPE svg>\n'+
     '<svg xmlns="http://www.w3.org/2000/svg"\n'+
     '  xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
     '  width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
    return s

  def plotend(sd) : return '</svg>\n'
  def move(sd,x,y) : return 'M'+str(x)+','+str(y)
  def line(sd,x,y) : return 'L'+str(x)+','+str(y)
      
  def labl(sd,s,x,y, size=14,fill='#000') :
    f=' font-family="Arial"'
    t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
    return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')

  def path(sd,w,s,lc='#000') : # path w=with lc=linecolor
    t='<path fill="none" stroke="'+lc+'" stroke-width="'+str(w)+'" '
    t+= ' stroke-linecap="round" ' 
    return t + 'd="\n'+s+'" />\n'

# end class svgdump

def aprx(x) : # 3-digit approximation
  y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
  return z/1000.0

def curve(sd,scale,xy, fill='#000', width=1) :
  (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2) 
  for i in range(n) : 
    x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
    s+= sd.move(x,y) if(i==0) else sd.line(x,y)
  return sd.path(width,s,lc=fill)

def label(sd,scale,txt,tx,ty) : # tx,ty plot coordinates not pixels
  (xs,fx,xmin, ys,fy,ymin) = scale 
  x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
  return sd.labl(txt,x,y)

###########

def shootlog(y,f,i,k) : # x=exp(y) version
  d=y[i]-y[i-1]; z=f[i]; j=i; p=1.0-0.5*d; q=1.0+0.5*d
  while (j<k) and (z>=0) and (z<10) : # stop if out of bounds
    t=exp(y[j])*z;  z=sqrt(t*t*t)
    u= d*d* z + 2*f[j]-q*f[j-1]
    f[j+1]= u/p; y[j+1]=y[j]+d; j+=1; z=f[j]
  return j,z  # print(str(f[j+1]))

def shootlin(x,f,i,k) : # i=1, f[0] and f[1] known, x[0]=0, x[1]=d.
  d=x[i]-x[i-1]; j=i; z=f[j]
  while (j<k) and (z>=0) and (z<10) :
    u= d*d* sqrt(z*z*z/x[j])
    f[j+1]= u+2*z-f[j-1]; x[j+1]=x[j]+d; j+=1; z=f[j]
  return j,z    # print(str(f[j+1]))

def logplot(fn,x,f,y,g) : # log x axis, y axis 0 to 1.
  txt=['Function y\"= sqrt(y^3/x), y(0)=1.',
    'xscale log 0.01...100','yscale 0...1']
  n=len(x); m=len(y); ok=(g[1]==f[n-1]); k=0; xy=[0.0]*(2*n+2*m)
  if not ok : print('logplot no join'); exit()
  for i in range(1,n) : xy[k]=log(x[i]);xy[k+1]=f[i]; k+=2
  for i in range(2,m) : xy[k]=y[i]; xy[k+1]=g[i]; k+=2
  xmin=log(0.01); xmax=log(100); ymin=0; ymax=1.0; xy=xy[:k]
  xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin ] # frame
  xg=[1,2,4,6,8]; fg=[0.01,0.1,1,10,100]; z=0; 
  lbl=['0.01','0.1','1','10','100']
  for i in range(len(fg)) :
    for k in range(len(xg)) :
      f=fg[i]; g=f*xg[k]; h=log(g)
      if (i==0)and(k==0) : xp=xmin; z=0
      elif (h<=xmax)and(z==0) : xyf += [xp,ymin,h,ymin,h,ymax];xp=h; z=1
      elif (h<=xmax)and(z==1) : xyf += [xp,ymax,h,ymax,h,ymin];xp=h; z=0
  xys=[xyf,xy]
  xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale user units, min in pixels
  ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#777','#000','#00a','#077','#0a0','#770','#a00']
  sd=svgdump(); buf=sd.plotbegin(800,600,1); i=0; dy=(ymax-ymin)/20.0
  for j in range(len(txt)) :
    buf += label(sd,scale, txt[j], (xmin+xmax)*0.5,ymax-dy-j*dy)
  for j in range(len(lbl)) :
    buf += label(sd,scale, lbl[j], log(fg[j]),ymax)
  for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

def shooter() : # dichotomic shots to solve equation
  # adjust f[1] to hit target
  # logscale y,g(y) after initial x,f(x) linear section 
  n=50; n1=n+1; x=[0]*n1; f=[0]*n1; x[0]=0;x[1]=0.01; f[0]=1; fhi=1; flo=0.5
  m=300; m1=m+1; y=[0]*m1; g=[0]*m1; logsection=True
  targ=1e-8; ff=0.5*(fhi+flo); f[1]=ff; zhi=10; zlo=-1; 
  i=0; busy=True;zbefore=1
  while busy and (i<100) :
    j,z= shootlin(x,f,1,n); print(str([f[1],j,z]))
    if logsection and (z>0)and(z<1) :
      g[0]=f[j-1]; g[1]=f[j]; y[0]=log(x[j-1]); y[1]=log(x[j])
      k,z= shootlog(y,g,1,m); xmax=exp(y[k]); print('log: '+str([xmax,k,z]))
    if z>targ : fhi=ff;zhi=z 
    else : flo=ff;zlo=z
    ff=0.5*(fhi+flo);
    if (zlo>0)and(zlo<zhi)and(zhi<1) : ff= flo+(targ-zlo)*(fhi-flo)/(zhi-zlo)
    f[1]=ff; i+=1; busy= (abs(z-zbefore)>1e-8); zbefore=z
  print('Done at i='+str(i))
  logplot('thomferm', x,f,y,g)
  atomradius('bla', x,f,y,g)

def atomradius(fn,x,f,y,g) : # solve x(a): u(x)-xu'(x)=a where a=0.01..0.1
  # u-xu'= u[i]-xx[i]*(u[i+1]-u[i-1])/(xx[i+1]-xx[i-1])
  n=len(x); m=len(y); k=0; xx=[0.0]*(n+m);u=[0.0]*(n+m); v=[0.0]*(n+m)
  for i in range(1,n) : xx[k]=x[i];u[k+1]=f[i]; k+=1
  for i in range(2,m) : xx[k]=exp(y[i]); u[k+1]=g[i]; k+=1
  vmax,jmax,vmin,jmin= -1e6,0,1e6,0
  for j in range(1,k-1) :
    v[j]= u[j]-xx[j]*(u[j+1]-u[j-1])/(xx[j+1]-xx[j-1])
    if v[j]>vmax: vmax,jmax= v[j],j
    if v[j]<vmin: vmin,jmin= v[j],j
  if False : #debug
    print('v='+str([v[1],v[50],v[100],v[k-2]]))
    print('vmax,xmax='+str([vmax,xx[jmax]]))
    print('vmin,xmin='+str([vmin,xx[jmin]]))
  print('Percent electron charge; Radius of ball')
  for p in range(1,21) : # p=percentage
    a=0.01*p ; j=jmax
    while v[j]>a : j+=1
    xa=xx[j-1];xb=xx[j]; va=v[j-1];vb=v[j]
    xp=xa+(a-va)*(xb-xa)/(vb-va); print('p,xp='+str([100-p,xp]))  

### following stuff not working

def geterror(f,y,err, n) :
  d=1.0/n
  for i in range(1,n) :
    x=y[i]; z=1.0-x; a=f[i]; f3=a*a*a
    f1=(f[i+1]-f[i-1])/(2.0*d); f2=(f[i+1]+f[i-1]-2*f[i])/(d*d)
    b=f2-2*z*f1; c=x*z*z*b*b; err[i]= f3-c

def shootstep(y,f,i,n) : # get f[i+1] for zero error.
  d=1.0/n; x=y[i]; z=1.0-x; a=f[i]; f3=a*a*a;
  # init bracket problem minimum find
  def err(i,ff) :
    f1=(ff-f[i-1])/(2.0*d); f2=(ff+f[i-1]-2*f[i])/(d*d)
    b=f2-2*z*f1; c=x*z*z*b*b; return f3-c
  # shoot is initialized with f[0] and f[1], the first slope
  def bugsweep(i,fh,eh) : # nearly degen pair of zeroes!
    df=1e-2; fmin=fh;emin=1e9
    for k in range(25) :
      ff=fh+(k-10)*df; ef=err(i,ff);
      fmin,emin = (ff,ef) if (abs(ef)<emin) else (fmin,emin)
      print('bug k,ff,ef='+str([k-10,ff,ef]))
    df=1e-3; fs=fmin; ok=False; fl=fmin;el=emin
    for k in range(25) :
      ff=fs+(k-10)*df; ef=err(i,ff)
      print('bug2 k,ff,ef='+str([k-10,ff,ef]))
      if (not ok) and ((ef*eh)<0) : ok=True; fl=ff;el=ef
    return fl,el
  #
  def descent(i,fh) :
    # search for bracket for f[i+1], fh<0 search a point >0
    eps=0.001; df= max(abs(eps*fh),eps); fl=fh-df; j=0; nf=0
    eh= err(i,fh); el=err(i,fl); found=(el*eh)<0; q=el/eh
    if (not found) and (q>1) : df=-df; fl=fh-df; el=err(i,fl)
    #print('i,eh,el,fh,fl='+str([i, eh,el, fh,fl]))
    while (j<100)and(not found) :
      q=el/eh
      #if q>1 : df=-0.5*df; fl -= df; nf+=1; print('flip j,eh,el='+str([j,eh,el]))
      if q>0.9 : df= 4*df; fl -=df
      elif q>0.5 : df= 2*df; fl -=df
      else : fl -= df
      el=err(i,fl); found=(el*eh)<0; j= 100 if (nf>4) else (j+1)
    if not found : fx,ex= bugsweep(i,fh,eh); 
    # print('i,j,eh,el,fh,fl='+str([i,j, eh,el, fh,fl]))
    return fl,el,j
  #
  fh=a; eh= err(i,fh); ok=True; k=0
  fl,el,j =descent(i,fh)
  while (k<20)and(ok) : # this part working, slow convergence
    if (eh*el)>=0 : ok=False;print('shoot exit i,k='+str(i)+','+str(k))
    ff= (fh*el-fl*eh) / (el-eh); ef=err(i,ff)
    if ef*eh>0 : fh=ff; eh=ef
    else : fl=ff; el=ef
    k +=1
  if ok : f[i+1]=ff; print('i,j,k, ff,ef='+str([i,j,k,ff,ef]))
  return ok

def smooth(x) :
  n=len(x); n1=n-1; y=[0]*n
  for i in range(1,n1) : y[i]=(2.0*x[i]+x[i-1]+x[i+1])/4.0
  return y

# def downerr(y,f,g,y,ef,eg, n) :
# tweak f[i],g[i] so that err at [i-1] doesnt worsen 

def explore(n) :
  d=1.0/n;n1=n+1; f=[0]*(n+1); g=[0]*n1; y=[0]*n1; f[0]=1; f[n]=0; 
  for i in range(n+1) : f[i]= 1.0-d*i; y[i]=d*i; g[i]=f[i] - 0.7*y[i]*f[i]
  print('start end', str([f[0],g[0],f[n],g[n]]))
  errf=[0]*n1; errg=[0]*n1; h=[0]*n1; errh=[0]*n1; h[0]=1
  geterror(f,y,errf, n); ef=smooth(errf); ef=errf
  geterror(g,y,errg, n); eg=smooth(errg); eg=errg
  for i in range(1,n) : h[i]=(f[i]*eg[i]-g[i]*ef[i]) / (eg[i]-ef[i])
  for i in range(1,n) : h[i]=(f[i]+g[i])/2.0
  geterror(h,y,errh,n) # secant fails, with eg,ef worse than average
  #for i in range(1,n) :
  #  print(str([i,errf[i],errg[i],errh[i]]))
  for i in range(1,n) : ok= shootstep(y,g,i,n)

# explore(20)
shooter()