Quantenmechanik/ H-Atom-Spektrum/ Pythonskript

Skript plplot.py Bearbeiten

'''
plplot.py  Legendre polys plot
l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
Pl(z):= (1/2^l l!) d^l (zz-1)^l
Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1

P0= 1
P1=z   P11=-x 
P2=(3zz-1)/2; P21= -3xz P22=3xx 
P3=(5zzz-3z)/2 P31=-3/2x(5zz-1) P32=15zxx; P33=-15xxx
P4=(35zzzz-30zz+3)/8
'''

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

def divide(p,q) : # improve integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  if p==0 : q=1
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return [s*p,q]

def polynscale(u,p,q) : # scale factor on polynomial
  n=len(u); z=[[]]*n
  for i in range(n) : z[i]= divide(p*u[i][0],q*u[i][1])
  return z

def mult(a,b) : # multiply rationals or integers, result rational
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return [pa*pb,qa*qb]

def add(a,b) : # add rational or int
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return divide(pa*qb+pb*qa,qa*qb)

def polynadd(u,v,a,b) : # rational-coeff polynomial added, with int? factors
  n=len(u); m=len(v); k=max(n,m); w=[[]]*k
  for i in range(k): 
    uu= mult(a,u[i]) if (i<n) else [0,1]
    vv= mult(b,v[i]) if (i<m) else [0,1]
    q=uu[1]*vv[1]; p=uu[0]*vv[1]+vv[0]*uu[1]; w[i]= divide(p,q)
  return w

def polynmul(u,v) : # multiply polynomials
  n=len(u); m=len(v)
  if (n==0)or(m==0) : w=[[0,1]]
  else :
    w=[[]]*(n+m-1)
    for i in range(n) :
      for k in range(m) :
        a=u[i][0];b=u[i][1]; c=v[k][0];d=v[k][1]; z= divide(a*c,b*d); j=i+k
        if w[j]==[] : w[j]=z # first use of j
        else : p=z[0]*w[j][1]+z[1]*w[j][0]; q=z[1]*w[j][1]; w[j]=divide(p,q)
  return w     

def polyndiff(u) : # derivative of polynomial
  n=len(u); w=[[]]*(n-1)
  for i in range(1,n) : w[i-1]= [i*u[i][0],u[i][1]]
  return w

def plzero(n) : # rational polys of z
  #  l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
  #  l Pl(z,x)= (2l-1) z P(l-1)(z,x) - (l-1)(zz+[1~=xx])P(l-2) 
  pl=[[]]*20; pl[0]= [[1,1]]; z= [[0,1],[1,1]]; un=[[1,1],[0,1],[1,1]]
  for i in range(1,n) :
    y= polynmul(pl[i-1],z); y= polynscale(y, 2*i-1,1)
    if (i-2)>= 0 : w= polynscale(pl[i-2],i-1,1); y=polynadd(y,w,1,-1)
    pl[i]= polynscale(y,1,i)
  return pl

def pleval(p,zmin,zmax,dz) : 
  nz= int((zmax-zmin)/(dz+0.0)) + 1; m=len(p); xy=[0]*(2*nz)
  for i in range(nz) :
    x=zmin+i*dz; y=0; zz=1
    for k in range(m) : y += zz* p[k][0]/(p[k][1]+0.0); zz *=x
    xy[2*i]=x; xy[2*i+1]=y
  return xy


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 miniplot(fn, xys, text) :
  colo=['#000','#000','#00a','#077','#0a0','#770','#a00']
  xs,fx,xmin= -1,  750/2.0 ,25 # xscale start in user units, min in pixels
  ys,fy,ymin= -1, -550/2.0 ,575 ; scale=(xs,fx,xmin, ys,fy,ymin)
  sd=svgdump(); buf=sd.plotbegin(800,600,1)
  buf += sd.labl(text,50,50); i=0
  for xy in xys : buf += curve(sd,scale,xy, fill=colo[i%7]); i +=1
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

def legendretest() :
  xyf=[1,-1, -1,-1, -1,1, 1,1, 1,-1, 0,-1, 0,1, 1,1, 1,0, -1,0]; xys=[xyf]  
  pl=plzero(6)
  for i in range(6) : xy=pleval(pl[i], -1,1, 0.01); xys += [xy]
  miniplot('plplot', xys, 'Legendre Pl_1 to Pl_5')

# radial funcs u(x)/x= x^l exp(-x) radpol(n,l,x) Laguerre-type polyns, n>l

def division(p,q) : #  integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  def div(p,q) : return int(p/q)
  pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return s*p,q

def radialpol(nmax=5) : # coeffs of radial hydrogen atom polyns
  print('Radialpolynome.'); pols=[]; # n=1,2,3,4, l=0..n-1, a0=1. 
  for n in range(1,nmax+1) :
    for l in range(n) :
      j=2*l+2; k=0; p=1; q=1; t=''; pol=[]
      while p != 0 :
        p,q= division(p,q); pol += [[p,q]]
        t += ' '+ (str(p) if(q==1) else (str(p)+'/'+str(q)))
        q *= (k+1)*(k+j); p *= (-2*n+2*k+j); k+=1 
      print('n='+str(n)+' l='+str(l)+' '+t)
      pols += [[n,l,pol]]
  return pols

def radiallist(nlist,mode=0) : # laguerre-type radial functions
  def pol(q,x) : # polynomial q
    y=0; zz=1; m=len(q)
    for k in range(m) : y += zz* q[k][0]/(q[k][1]+0.0); zz *=x
    return y
  testrun= (mode==0)
  p=radialpol(); m=len(p);
  dx=0.05; xmax=15.0;xmin=0.0; ymin=-1;ymax=1; xys=[]; minmax=[]; k=0
  for i in range(m) :
    n=p[i][0]; l=p[i][1]; q=p[i][2]; 
    if n==nlist : # scan x^l expm(x) pol(q,x), 200 points xy[l]
      x=0; y1=0;y=0; j=0; ymi=0;yma=0
      while (x<1)or(abs(y)>1e-3)or(y1>1e-3) :
        z= exp(-x)*pol(q,x)
        for k in range(l) : z *=x
        y1=abs(y); y=z; x +=dx; j+=1; ymi=min(ymi,y); yma=max(yma,y)
      #print('n,l,j,min,max='+str([n,l,j,ymi,yma])); minmax += [(ymi,yma)]
      if not testrun :
        x=0; xy=[]; a,b= ymi,yma; b= (-a) if((-a)>b) else b; b=1.01*b
        while x <= xmax :
          z= exp(-x)*pol(q,x)
          for k in range(l) : z *=x
          xy += [x,z/b]; x += dx
        xys += [xy]
      k+=1 # after defining xys[k] and/or minmax[k] 
  xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin,
        xmin,0, xmax,0] # frame    
  xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale start user units, min in pixels
  ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#000','#000','#00a','#077','#0a0','#770','#a00']
  fn='radialf'; text='Radial n=5 l=0,1,2,3,4'; xys= [xyf] + xys
  if not testrun : # strange plot, case l=0 tiny, l=4 huge.
    sd=svgdump(); buf=sd.plotbegin(800,600,1)
    buf += sd.labl(text,50,500); i=0
    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()

legendretest()
radiallist(5,1)    

Skript pltable.py Bearbeiten


# Legendre polynomials, rational-number arith

def div(p,q) : return int(p/q)

def getprimes(n) : # want n prime numbers
  p=[2]; i=1; plast=2
  while i<=n :
    q=plast+1; found=False
    while not found:
      j=0
      while ((q%p[j])!=0) and ((p[j]*p[j])<=q) :  j +=1
      if (q%p[j] != 0) : p += [q]; plast=q; found=True; i+=1
      else : q+=1
  return p

def greatestdiv(x) : # greatest divisor, skip 0s in list
  n=len(x); y=[0]*n; w=[0]*n; h=0; ymin=-1; gcd=1
  pn=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]; ln= len(pn) # small primes
  for i in range(n) :
    z=abs(x[i])
    if z>0 : ymin= z if(ymin<0) else min(ymin,z); w[h]=z; h+=1
  w=w[:h]; ip=0
  while (ip<ln)and(pn[ip]<=ymin) :
    z=pn[ip]; ok=True
    while ok :
      for i in range(h) : ok= ok and ((w[i]%z)==0)
      if ok : 
        for i in range(h) : w[i]= div(w[i],z)
        ymin= div(ymin,z); gcd *= z
    ip +=1 
  return gcd    

def smallestmult(x) : # smallest multiple of nonzero numbers
  n=len(x); y=abs(x[0])
  for i in range(1,n) : 
    z=abs(x[i]); y = y if (z==0) else div(y*z,greatestdiv([y,z]))
  return y

def divide(p,q) : # improve integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  if p==0 : q=1
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return [s*p,q]

def intsolve(ai, bi) :  # solve integer equation, return fractions as pairs
  # ai is p-by-p matrix, bi is p-vector.  solve Ax=b 
  def fracadd(x,y,f) : # x + f*y . x,y fractions, f integer.
    a,b= tuple(x); c,d= tuple(y); c *=f; return divide(a*d+c*b,b*d)
  p=len(bi); xg=[[]]*p; bg=[0]*p; ok= True; ag=[[]]*p; j=0 
  for i in range(p) : # local copy
    ag[i]=[0]*p; bg[i]= int(bi[i])
    for k in range(p) : ag[i][k] = int(ai[i][k])
  while ok and (j<p) : # kill lower triangle columns, row by row
    i=j # get pivot. swap first non-zero ag[j][j] row up
    while (ag[i][j]==0) and (i<p) : i +=1
    ok= (i<p)
    if ok and (i>j) : 
      for h in range(p) : z=ag[i][h]; ag[i][h]=ag[j][h]; ag[j][h]=z
      z=bg[i]; bg[i]=bg[j]; bg[j]=z
    for i in range(j+1,p) :
      rp= ag[i][j]; rq= ag[j][j] # line i = rq*(line i) - rp* (line j) 
      bg[i]= rq*bg[i]-rp*bg[j]; ag[i][j]= 0; ok= ok and (rq!=0)  
      for l in range(j+1,p) : ag[i][l]= rq*ag[i][l]-rp*ag[j][l] 
    j +=1
  if ok : #  matrix ag is upper triangular 
    i=p-1
    while i>=0 : # make vector xg: (ag) * (xg) = (bg) 
      ax= [0,1]
      for j in range(i+1,p) : ax= fracadd(ax,xg[j],ag[i][j]) 
      z= fracadd([bg[i],1],ax,-1); xg[i]= divide(z[0],z[1]*ag[i][i]); i -=1 
  return ok, xg # junk if ok=False

def polynscale(u,p,q) : # scale factor on polynomial
  n=len(u); z=[[]]*n
  for i in range(n) : z[i]= divide(p*u[i][0],q*u[i][1])
  return z

def mult(a,b) : # multiply rationals or integers, result rational
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return [pa*pb,qa*qb]

def add(a,b) : # add rational or int
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return divide(pa*qb+pb*qa,qa*qb)

def polynadd(u,v,a,b) : # rational-coeff polynomial added, with int? factors
  n=len(u); m=len(v); k=max(n,m); w=[[]]*k
  for i in range(k): 
    uu= mult(a,u[i]) if (i<n) else [0,1]
    vv= mult(b,v[i]) if (i<m) else [0,1]
    q=uu[1]*vv[1]; p=uu[0]*vv[1]+vv[0]*uu[1]; w[i]= divide(p,q)
  return w

def polynmul(u,v) : # multiply polynomials
  n=len(u); m=len(v); w=[[]]*(n+m-1)
  for i in range(n) :
    for k in range(m) :
      a=u[i][0];b=u[i][1]; c=v[k][0];d=v[k][1]; z= divide(a*c,b*d); j=i+k
      if w[j]==[] : w[j]=z # first use of j
      else : p=z[0]*w[j][1]+z[1]*w[j][0]; q=z[1]*w[j][1]; w[j]=divide(p,q)
  return w     

def polyndiff(u) : # derivative of polynomial
  n=len(u); w=[[]]*(n-1)
  for i in range(1,n) : w[i-1]= [i*u[i][0],u[i][1]]
  return w

def wholenb(x) : # x is a list rationals, scale to get integers only. bad.
  n=len(x); q=1; y=[0]*n
  for i in range(n) : r=x[i][1]; z=divide(q,r); q = z[0]*r
  for i in range(n) : z=divide(q*x[i][0],x[i][1]); y[i]= z[0]
  if y[0]<0 :
    for i in range(n) : y[i]= -y[i]
  return y

def polynpurge(u) : # omit leading zero terms and scale to integers
  # return new polyn and scalefactor applied
  n=len(u); j=n-1; k=0 
  while (u[j][0]==0)and (j>=0) : j -=1
  while (u[k][0]==0)and (k<n) : k+=1
  j +=1; v= u[:j]; p=[0]*j; q=[0]*j; sign= (-1) if (u[k][0]<0) else 1
  for i in range(j) : p[i]= v[i][0]; q[i]=v[i][1]
  z=smallestmult(q) 
  for i in range(j) : r= divide(p[i]*z,q[i]); p[i]=r[0]
  y=greatestdiv(p); 
  for i in range(j) : r=divide(p[i],y); v[i]=[sign*r[0],1]  
  return v, divide(z,y)

def present(s,q,mode,noisy) :
  # mode 0 no x, mode 1 complementary x, mode -n const x^n, mode>0 l=mode 
  n=len(q); dmax=1; j= (n-1) if (mode==1) else 0; s+='('
  if mode>0 : j=mode # new
  v= '' if (mode >=0) else ('x'+str(-mode))
  for i in range(n) : dmax= q[i][1] if (q[i][1]>dmax) else dmax
  for i in range(n) :
    a=q[i][0]; b= q[i][1]; bad= (dmax % b) != 0; sig=''
    if bad : print('present bug exit'); exit()
    t=('x'+str(j)) if (j>0) else ''; u=('z'+str(i)) if (i>0) else ''
    if a!=0 : cf= div(dmax,b)*a; sig='+' if (cf>0) else ''
    if a!=0 : s+= ' '+sig+str(cf)+v+t+u  
    j -=1
  s+=')/'+str(dmax) 
  if noisy : print(s)
  return (s+'\n') 

def plhomogenize(q,l,i) : # degree l, implicit x-power i. 
  # multiply with (1+zz) truncated poly up to k, again up to k-2 etc.
  # powers x are implied, the complementary ones to z
  # inverse plreduce on sphere: drive out all xx as (1-zz) leaving only x.
  n=len(q); f=[[1,1],[0,1],[1,1]] # factor polyn 1+zz to apply
  t= []; ok=True; dg=[0]*n # t target polyn 
  for k in range(n) :
    if q[k][0] != 0 : j=l-(k+i); dg[k]=div(j,2); ok= ok and ((j%2)==0) 
  if not ok : print('plhomogenize bad input. exit.'); exit()
  #print('dg='+str(dg)) # dbg
  for k in range(n) :
    u=[[0,1]]*(k+1); u[k]=[]+q[k]; h=dg[k]
    # if h>0 : multiply coeff with f^h, then add to t
    for j in range(h) : u=polynmul(u,f)
    t= polynadd(t,u,1,1)
  t,f= polynpurge(t); return t # divide out common factors 

def plmore(pl,l,noisy) :
  #Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1
  # option: homogenize with products (xx+zz) when deficient
  ql=pl;pwr=1;sign=-1; hm='';pi=''; rl=[[]]*(l+1);rl[0]=[]; qlist=[[l,0,[]+pl]]
  for i in range(1,l+1) :
    ql= polyndiff(ql); pwr += 1; ql=polynscale(ql,sign,1); qlist+=[[l,i,ql]]
    pi += present('Plm['+str(l)+str(i)+']=', ql, -i, noisy)
    # power i of x included, z^k must get (l-(k+i))/2 factors (1+zz). 
    rl[i]= plhomogenize(ql,l,i)
    hm += present('Hlm['+str(l)+str(i)+']=', rl[i], l, noisy)
  return hm,pi,rl,qlist
 
def plzero(n,noisy) : # rational polys of z
  #  l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
  #  l Pl(z,x)= (2l-1) z P(l-1)(z,x) - (l-1)(zz+[1~=xx])P(l-2) 
  # returns hom,hpl: homog versions in (x,z); pl,pli: official Plm
  # have either odd-only or even-only powers. let p be the maxpower
  # multiply term p-2 with (xx+zz), (p-2k) with (xx+zz)^k.
  # store in same array type, powers of x are implicit.
  pl=[[]]*20; pl[0]= [[1,1]]; z= [[0,1],[1,1]]; un=[[1,1],[0,1],[1,1]]
  ql=[[]]*20; ql[0]= [[1,1]]; pli='Plm[0,0]=(1)\n'; pllist=[]
  hom=''; hpl=[[[1,1]]] # homogeneous variant, stringlist and poly list
  for i in range(1,n) :
    y= polynmul(pl[i-1],z); y= polynscale(y, 2*i-1,1)
    if (i-2)>= 0 : w= polynscale(pl[i-2],i-1,1); y=polynadd(y,w,1,-1)
    pl[i]= polynscale(y,1,i)
    pli += present('Plm['+str(i)+'0]=',pl[i], 0, noisy)
    y= polynmul(ql[i-1],z); y= polynscale(y, 2*i-1,1) # homog version
    if (i-2)>= 0 : # wrong for i>2
      w= polynscale(ql[i-2],i-1,1); w= polynmul(w,un); y=polynadd(y,w,1,-1)
    ql[i]= polynscale(y,1,i); qq,fq = polynpurge(ql[i]) # must keep ql[i] as-is
    hom += present('Hlm['+str(i)+'0]=',qq, i,noisy)
    hm,pi,rl,qlist = plmore(pl[i],i, noisy)
    pllist += qlist; rl[0]=qq;  hpl += rl
    hom+=hm;pli+=pi # expand powers x^(2n) as (1-zz)^n
  return pl,hom,pli,hpl,pllist

def listpos(lst,i,j,k) :
  n=len(lst); ix=-1
  for h in range(n) : ix=h if([i,j,k]==lst[h]) else ix
  if ix<0 : lst += [[i,j,k]]; ix=n
  return lst,ix  

def combine(l,m,base) :
  n=len(base); cf=[[]]*(2*n); cfs=[]; icf=0
  for h in range(n) :
    i,j,k=tuple(base[h]); i1,j1,k1= i-1,j-1,k; i2,j2,k2= i,j,k-2; cfok=True
    ok1= (i1<0)or(j1<0); ok2= (k2<0) # vanishing terms
    if not ok1 : cf[icf]= [4*i*j,i1,j1,k1,h]; icf+=1 
    if not ok2 : cf[icf]= [k*(k-1),i2,j2,k2,h]; icf +=1 
  if n==1 : # check that delta deriv terms vanish, ie have a negative ix
    i,j,k=tuple(base[0]); cfs=[1]
    i1,j1,k1= i-1,j-1,k; i2,j2,k2= i,j,k-2
    ok= ((i1<0)or(j1<0)) and (k2<0)
    if not ok : print('combine bug: '+str([i1,j1,k1])+str([i2,j2,k2]))
  if (n==2)and False : # must have 1 match and 2 negatives, old code
    for h in range(1,4): cfok= cfok and (cf[0][h]==cf[1][h])
    # in general extract all cf with same 123 triplets, return matrix
    cfs=[cf[1][0],-cf[0][0]] # winning zero lincomb 
    print('combine cfs:'+str(cfok)+' '+str(cf)+' '+str(cfs)) 
  if n>=2 : # make equation system, index different ijk with position p
    mx=[[]]*n; lst=[]; b=[0]*n # mx has 2 nontrivial lines, common center elem?
    for h in range(n) : mx[h]=[0]*n
    for h in range(n) : mx[n-1][h]=1; b[n-1]=1 # dummy row
    for h in range(icf) :
      i,j,k= tuple(cf[h][1:4]); 
      lst,p= listpos(lst,i,j,k); q=cf[h][4] # q= basis ix
      # make n*n matrix so that coeff sums for same p become zero
      # add a line with maxvals of cols? or sum=1?
      mx[p][q]= cf[h][0]
    # print('matrix '+str(mx)); # showmatrix(m) 
    ok,xg=intsolve(mx,b); # boost to smallest ints and return as cfs 
    cfs= wholenb(xg)
    if False : print('intsolve ok='+str(ok)+' '+str(xg)+' cfs='+str(cfs))
  return cfs

def listterms(lmax,mode,noisy) :
  # loop over (l,m),  monomials of Ylm to polynomials
  # in terms of Plm exp(im f), Plm(x,z) is obtained by x^(i+j) z^k
  # when x= sin(teta); z= cos(teta) this written in c s smbols?
  # note Pl(-m)=Plm: because i-j=0 iff j-i=-m.
  full=(mode==1); data=[]; text=''; dm=-1
  for l in range(lmax+1) :
    m=l; lastm= (-l) if full else 0
    if mode==2 : m=0; dm=1; lastm=l # growing order
    while (dm*(lastm-m)) >= 0 :
      i=l; s='l='+str(l)+' m='+str(m); base=[]; t=''+s; u=''+s; block=[l,m]
      while i>=0 :
        j=i-m; k=l-i-j
        if (j>=0)and(k>=0) : s += ' '+str([i,j,k]); base += [[i,j,k]]
        i-=1       
      cfs= combine(l,m,base); lb=len(base)
      for h in range(lb) : # new neutral: cc becomes x, ss becomes z
        z=str(cfs[h]); z= ('+'+z) if (cfs[h]>0) else z; i,j,k=tuple(base[h])
        g=i+j; cc='' if (g==0) else ('x' if (g==1) else 'x^'+str(g)) 
        ss='' if (k==0) else ('z' if (k==1) else 'z^'+str(k)); void=(cc+ss)=='' 
        if not void : z= '+' if (cfs[h]==1) else ('-' if (cfs[h]==-1) else z)
        t += ' '+z+str(base[h]); u+= ' '+z+cc+ss; block+=[base[h]+[cfs[h]]]  
      data += [block]
      if noisy : print(u)
      m += dm; text += u +'\n'
  return data,text

'''
  I(m,n)= integ 0 to pi of sin^m cos^n [ n even ]
  - m += 2 implies  I(m+2,n) = (m+1)/(m+n+2) *I(m,n)
  - n += 2 implies  I(m,n+2) = (n+1)/(m+n+2) *I(m,n)
  - start functions : 1, sin, cos, sincos: 0,0 1,0 0,1, 1,1 
  0,0: pi. 1,0: 2 0,1 : 0  1,1 : 0
'''
def isico0pi(m,n) : # integral 0 to pi of sin^m(x) cos^n(x)
  if (n%2 != 0) : return [0,1] # ugly trick: negative val encode pi factor
  nn=0; mm= 0 if (m%2 == 0) else 1; p=-1 if (mm==0) else 2; q=1 
  while nn<n : p *= (nn+1); q *= (mm+nn+2); nn += 2
  p,q= tuple(divide(p,q))
  while mm<m : p *= (mm+1); q *= (mm+nn+2); mm += 2
  return divide(p,q)   

def csproduct(n,a,b,c, m,d,e,f) : # product of 2-var polynomials p(x,z)
  # length n and m. a,d= lists of coeffs, b,e powers of x, c,f powers of z.
  # returns list of coeffs and list of power pairs they target
  # isico0pi uses these power pairs, results to be summed up then
  # variant: a,d are pairs=rationals.
  def ratpr(a,b) : return [a[0]*b[0],a[1]*b[1]]
  def ratsum(a,b) :
    if a==[] : a=[0,1]
    p=a[0]*b[1]+a[1]*b[0];q=a[1]*b[1]; return divide(p,q)
  def getindex(ix,p,q) : # where to store x^p y^q
    i=0; n=len(ix); found=False
    while (not found) and (i<n) and (ix[i]!=[]) :
      found= (ix[i][0]==p)and(ix[i][1]==q); i= (i+1) if (not found) else i
    if not found : ix[i]=[p,q] # new exponent combination
    return i 
  t=str(type(a[0])); ratnum= (t.find('list')>0)
  x= ([[]]*(n+m)) if ratnum else ([0]*(n+m)) 
  ix=[[]]*(n+m); jmax=0
  for i in range(n) :
    for k in range(m) :
      yy=b[i]+e[k]; zz=c[i]+f[k]; j=getindex(ix,yy,zz) # side effect on ix
      if ratnum : x[j]= ratsum(x[j],ratpr(a[i],d[k]))
      else : x[j] += a[i]*d[k]
      if j>jmax : jmax=j
  return x[:(jmax+1)],ix[:(jmax+1)]

def intprod(pa,ia,la, pb,ib,lb) : # integral of product of 2 plm(x,z) polyn
  # pa(x,z)= sum(k,pa[k]*x^ia*z^k, pb same. nonhomog, const x power=i+j
  # status ok for plzero
  na=len(pa); cfa=[0]*na; pxa=[0]*na; pza=[0]*na
  for i in range(na) : cfa[i]=pa[i]; pxa[i]=ia; pza[i]=i
  nb=len(pb); cfb=[0]*nb; pxb=[0]*nb; pzb=[0]*nb
  for i in range(nb) : cfb[i]=pb[i]; pxb[i]=ib; pzb[i]=i
  cf,pw= csproduct(na,cfa,pxa,pza, nb,cfb,pxb,pzb); nt=len(cf); sum=[]
  # print('intprod cf,pw='+str([cf,pw]))
  for j in range(nt) :
    pw[j][0] +=1 # sin factor integral measure
    z= isico0pi(pw[j][0],pw[j][1]); a=z[0]; b=z[1]
    z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
    sum= polynadd(sum,z, 1, cf[j])
  # print('pa,sum='+str([pa,sum]))
  return sum

def iproducts(dat,all,noisy) :
  # got listterms data for Plm, eval int(Plm^2 sin(a)da)
  # check, int(Pl1m1 Pl1ml2m sin(a) da) must vanish if l1 != l2.
  # scalar prods either have a factor pi or not, in all terms.
  mode=2 # mode2: cs are sinuses, sn cosinuses. Plm(x,z): z~cos, x~sin.
  n=len(dat); csh=[0]*32; snh=[0]*32; cfh=[0]*32; bugs=0
  csk=[0]*32; snk=[0]*32; cfk=[0]*32; tbl='' # tbl result table of squares
  for h in range(n) :
    x=dat[h]; l=x[0];m=x[1]; jh=len(x)-2; kmax= n if all else (h+1)
    for i in range(jh) : # cf=coeff, cs=cosinus and sn=sinus powers of Plm
      y=x[2+i]; cfh[i]=y[3]; csh[i]=y[0]+y[1]; snh[i]=y[2]
    for k in range(h,kmax) : # (h,h+1) or (h,n)  
      x=dat[k]; lk=x[0];mk=x[1]; jk=len(x)-2; sum=[]
      for i in range(jk) :
        y=x[2+i]; cfk[i]=y[3]; csk[i]=y[0]+y[1]; snk[i]=y[2]
      cf2,pw2= csproduct(jh,cfh,csh,snh,jk,cfk,csk,snk); nt=len(cf2) # nt terms
      for j in range(nt) : pw2[j][0] += 1 # increment sin-power of integrand
      for j in range(nt) : # arg order: sin-power, then cos-power
        if mode==1 : # bad
          z=isincos(pw2[j][1],pw2[j][0],0); sum= polynadd(sum,z, 1, cf2[j])
        elif mode==2 : # better integrals, sin/cos swapped, 0 to pi. ok!
          z= isico0pi(pw2[j][0],pw2[j][1]); a=z[0]; b=z[1]
          z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
          sum= polynadd(sum,z, 1, cf2[j])
      #print('l,m='+str([l,m])+' cf2='+str(cf2)
      # +' pw2='+str(pw2)+' sum='+str(sum))
      if noisy : print('l,m,lk,mk,nt='+str([l,m,lk,mk,nt])+' sum='+str(sum))
      z= sum[0][0] if (len(sum)==1) else (abs(sum[0][0])+abs(sum[1][0]))
      if (k==h) and (z==0) : print('BUG1'); bugs+=1 
      if k==h : tbl += 'l='+str(l)+' m='+str(m)+' '+str(sum)+'\n'
      if (l != lk) and (m==mk) and (z != 0) : print('BUG2'); bugs +=1
  if noisy : print('Bugs:'+str(bugs))
  return tbl
    
def getdecimal(s,n,trig) :
  d='0123456789'; i=0; m=len(s); p=[0]*n; t=''
  if trig != '' : i=s.find(trig)
  for k in range(n) :
    while (i<m) and (d.find(s[i])<0) : i+=1
    while (i<m) and (d.find(s[i])>=0) : t+=s[i]; i+=1
    p[k]= int(t); t=''
  return p

def stripblanks(s) :
  n=len(s); t=''
  for i in range(n) : t += '' if (s[i]==' ') else s[i]
  return t

def plmtable(txt) : # latex-style table
  u= txt.split('\n'); lu=len(u); v0='---'; v=[[]]*lu
  if u[lu-1]=='' : lu-=1; u= u[:lu]
  for i in range(lu) : v[i]= u[i].split(' ')
  t='\\begin{array}{|l|l|l|l|l|}\n'
  t+= ('\\ell & m & P_{\\ell m}(x,z) & P^h_{\\ell m}(x,z) homogen & '+
    '\\int_0^{\\pi} d\\theta\\; \\sin(\\theta)\\; (P^h_{\\ell m}'+
    '(\\sin(\\theta),\\cos(\\theta))^2 \\\\\n')
  for i in range(lu) :
    if (v[i][0] != v0) : t+= '\\hline\n'
    t+=v[i][0]+' & '+v[i][1]+' & '+v[i][2]+' & '+v[i][3]+' & '+v[i][4]+'\\\\\n'
    v0=v[i][0]
  t+= '\\end{array}\n'
  print(t); return(t)

def integlm(l,m) : # theoretical integral Plm^2, defining non-homog versions
  p=2; q=2*l+1; j=l-m+1
  while j<=(l+m) : p *= j; j +=1
  # print('integlm l,m,p,q'+str([l,m,p,q]))
  return divide(p,q)

def pleval(p,l,m) : # values at z=+-1
  if m>0 : return [] # any power of x=0 
  zv=[1,-1]; nz=len(zv); m=len(p); yv=[[]]*nz
  for i in range(nz) :
    z=zv[i]; y=0; zz=1
    for k in range(m) : y=add(y,mult(zz,p[k])); zz *=z
    yv[i]=y
  return yv

def testplm1(lmax=5) : # Plm build and table listing
  noisy=False
  plz,r,ri,hpl,pllist= plzero(lmax+1, noisy)
  if noisy : print(r);print('') # recursion
  d,txt = listterms(lmax,2, noisy) # lin-equation method
  tbl=iproducts(d,True,noisy); # integrals, odd trigon powers no pi factor.
  tx2=txt.split('\n'); tb2= tbl.split('\n'); n=len(tx2); m=len(tb2); dat=''
  for i in range(n) : 
    u=tx2[i]; s=tb2[i]
    if len(s) > 2 :
      j=u.find('+'); v=u[j+1:]; w=u[:j]; v=stripblanks(v)
      p= getdecimal(s,2,'[['); tb2[i]=(w+v+' '+str(p[0])+'/'+str(p[1]))
  ri2=ri.split('\n'); m=len(ri2)
  for i in range(m):
    u=ri2[i]
    if len(u)>0 : 
      a=u.find('('); v=stripblanks(u[a:]); b=v.find(')'); w='';c='?'
      if v[b+1:]=='/1' : v=v[:b+1]
      if v[1]=='+' : v=v[0]+v[2:]
      for k in range(len(v)) :
        block= (c=='x')or(c=='z');
        if block : w+= c if(v[k]=='1') else (c+'^'+ v[k])
        c=v[k]; w+= '' if ((c=='x')or(c=='z')or block) else c
      ri2[i]=w; # print(tb2[i]+' '+ri2[i]) # remix to final string 
  for i in range(m) :
    if len(ri2[i])>0 :
      q= tb2[i].split(' '); # print('['+q[0]+'|'+q[1]+']')
      x= getdecimal(q[0],1,''); y= getdecimal(q[1],1,'')
      s=str(x[0])+' '+str(y[0])+' '+ri2[i]+' '+q[2]+' '+q[3]
      dat= s if (dat=='') else (dat+'\n'+s)
  plmtable(dat)

testplm1()