' simpvai.bas - jack@ord.ca - 4 May 05 - Liberty Basic v4.01
' the algorithm optimizes the fit of P-A data a varying angle of incidence
'   to a model consisting of one or two layers on a known substrate
' NOTE: THE ALGORITHM DOES ARITHMETIC, BUT MAKES NO CLAIM ABOUT
'       THE SIGNIFICANCE OF THE PARAMETER VALUES IT DETERMINES

' Initialize Window
  nomainwin
  textbox#1.t1, 125, 3, 90, 20:   textbox#1.t4, 285, 3, 90, 20
  textbox#1.t2, 125, 33, 90, 20:  textbox#1.t5, 285, 33, 90, 20
  textbox#1.t3, 125, 63, 90, 20:  textbox#1.t6, 285, 63, 90, 20
  checkbox#1.c1, "Nin", [waitHere], [waitHere], 75, 3, 50, 20
  checkbox#1.c2, "Kin", [waitHere], [waitHere], 75, 33, 50, 20
  checkbox#1.c3, "Lin", [waitHere], [waitHere], 75, 63, 50, 20
  checkbox#1.c4, "Nout", [waitHere], [waitHere], 230, 3, 50, 20
  checkbox#1.c5, "Kout", [waitHere], [waitHere], 230, 33, 50, 20
  checkbox#1.c6, "Lout", [waitHere], [waitHere], 230, 63, 50, 20
  button#1, "Data1", [b1], UL, 3, 3, 50, 20
  button#1, "Data2", [b2], UL, 3, 28, 50, 20
  button#1, "Tplot", [b3], UL, 3, 53, 50, 20
  button#1, "Smplx", [b4], UL, 3, 78, 50, 20
  button#1, "Print", [b5], UL, 3, 103, 50, 20
  graphicbox#1.g1, 390, 3, 200, 80
  WindowWidth=610     ' pixel scale 0-600
  WindowHeight=459    ' pixel scale 0-430
  UpperLeftX=10: UpperLeftY=10
  open "Simplex algorithm for data at varying angle-of incidence" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim ip(8), db(13), q(11, 8), da(20, 4)
  dim c(4, 2), d(4, 2)
  dim b(2), u(2), v(2), j(4), k(4), m(4), h(4), r(4), s(4), e(4), f(4)
  pi=4*atn(1): kk=0: w=0: sn=7: f$="##.######": g$="######.##"
  #1 "down ; place 90 90 ; box 580 390"

[waitHere]
  wait

[sim]
  #1.g1 "cls ; down ; color black"
  for i=1 to 3
    x1=50*i: #1.g1 "line "; x1; " 0 "; x1; " 80"
  next i
  #1.g1 "color red"
  i=0: jlof=0
  tt=time$("ms" )
  dni=.02: dnf=.0000001
  if w=0 then zz$="Weight Off" else zz$="Weight On"
  #1.c1, "value?": input#1.c1, z$: if z$="set" then i=i+1: ip(i)=1
  #1.c2, "value?": input#1.c2, z$: if z$="set" then i=i+1: ip(i)=2
  #1.c3, "value?": input#1.c3, z$: if z$="set" then i=i+1: ip(i)=3
  if nlay=2 then
    #1.c4, "value?": input#1.c4, z$: if z$="set" then i=i+1: ip(i)=4
    #1.c5, "value?": input#1.c5, z$: if z$="set" then i=i+1: ip(i)=5
    #1.c6, "value?": input#1.c6, z$: if z$="set" then i=i+1: ip(i)=6
  end if
  nvr=i: nvx=i+1: ip(nvx)=sn
  if nvr=0 then
    df=0: iv=1: gosub [fit]
  else
    for i=2 to nvx+3
      for j=1 to 4
        q(i, j)=q(1, j)
      next j
    next i
    for i=1 to nvr
      if ip(i)=3 or ip(i)=6 then
        q(i, ip(i))=q(i, ip(i))+1
      else
        q(i, ip(i))=q(i, ip(i))+dni
      end if
    next i
    for iv=1 to nvx
      gosub [fit]
    next iv
    gosub [sort]
    y1=0
  end if
  while df>dnf
    x1=int(15*log(df/dnf)): y1=y1+1: if x1>190 then x1=190
    if y1=1 then z$="place " else z$="goto "
    #1.g1 z$; x1; " "; y1
    if y1=79 then y1=0
    xflag=0
    for i=1 to nvr
      ave=0
      for j=1 to nvr
        ave=ave+q(j, ip(i))
      next j
      q(nvx+1, ip(i))=ave/nvr
    next i
    if (jlof=1) and (q(nvx+2, sn)=q(nvx, sn)) then
      jlof=0: gosub [cns]
    else
      jlof=0: iv=nvx+2
      for i=1 to nvr
        q(iv, ip(i))=2*q(nvx+1, ip(i))-q(nvx, ip(i))
      next i
      gosub [fit]
      if q(iv, sn)>q(nvx, sn) then
        gosub [cns]
      else
        if q(iv, sn)<q(1, sn) then
          iv=nvx+3: nexp=2
          while xflag=0
            for i=1 to nvr
              q(iv, ip(i))=(nexp+1)*q(nvx+1, ip(i))-nexp*q(nvx, ip(i))
            next i
            gosub [fit]
            if q(iv, sn)<q(iv-1, sn) then
              for i=1 to nvx
                q(iv-1, ip(i))=q(iv, ip(i))
              next i
            else
              xflag=1
            end if
          wend
          for i=1 to nvx
            q(nvx, ip(i))=q(nvx+2, ip(i))
          next i
        else
          jlof=1
          for i=1 to nvx
            q(nvx, ip(i))=q(nvx+2, ip(i))
          next i
        end if
      end if
    end if
    gosub [sort]
  wend
  #1 "color blue"
  #1 "place 120 410": #1 "\rms dev   "; using(f$, q(1, sn))
  tt=time$("ms")-tt
  #1 "place 400 410": #1 "\exec time "; tt/1000; " sec"
  db(6)=da(ndata, 2): gosub [tplot]
  #1.t1, using(f$, q(1, 1))
  #1.t2, using(f$, q(1, 2))
  #1.t3, using(g$, q(1, 3))
  if nlay=2 then
    #1.t4, using(f$, q(1, 4))
    #1.t5, using(f$, q(1, 5))
    #1.t6, using(g$, q(1, 6))
  end if
goto [waitHere]

[sort]
  s=1
  while s>0
    s=0
    for i=1 to nvr
      if q(i, sn)>q(i+1, sn) then
        s=s+1
        for j=1 to nvx
          temp=q(i, ip(j)): q(i, ip(j))=q(i+1, ip(j)): q(i+1, ip(j))=temp
        next j
      end if
    next i
  wend
  for i=1 to nvr
    max=q(1, ip(i)): min=max
    for j=2 to nvx
      if (q(j, ip(i))>max) then max=q(j, ip(i))
      if (q(j, ip(i))<min) then min=q(j, ip(i))
    next j
    if i=1 then
      df=max-min
    else
      if (max-min)>df then df=max-min
    end if
  next i
return

[cns]
  iv=nvx+2
  for i=1 to nvr
    q(iv, ip(i))=(q(nvx+1, ip(i))+q(nvx, ip(i)))/2
  next i
  gosub [fit]
  if q(iv, sn)<q(nvx, sn) then
    for i=1 to nvx
      q(nvx, ip(i))=q(iv, ip(i))
    next i
  else
    for i=1 to nvr
      for j=2 to nvx
        q(j, ip(i))=(q(1, ip(i))+q(j, ip(i))) / 2
      next j
    next i
    for iv=2 to nvx
      gosub [fit]
    next iv
  end if
return

[dev]
  gosub [pa]
  if (pp-da(j, 0))>90 then pp=pp-180
  if (da(j, 0)-pp)>90 then pp=pp+180
  pd=da(j, 0)-pp
  if w>0 then pd=pd*sin(aa*pi/90)
  ad=da(j, 1)-aa: ds=pd*pd+ad*ad
return

[fit]
  fs=0
  for ij=1 to ndata
    theta=da(ij, 2): gosub [pa]
    if pp-da(ij, 0)>90 then pp=pp-180
    if da(ij, 0)-pp>90 then pp=pp+180
    ddp=da(ij, 0)-pp
    if w>0 then ddp=ddp*sin(aa*pi/90)
    dda=da(ij, 1)-aa: ds=ddp*ddp+dda*dda
    da(ij, 3)=sqr(ds): fs=fs+ds
  next ij
  q(iv, sn)=sqr(fs/ndata)
return

[b1]                                                   ' Data1 parameters
  nlay=1
  nsub=1.74969: ksub=0: namb=1.00027
  db(1)=1.46: db(2)=0: db(3)=1148
  db(7)=45: db(8)=180: db(9)=0: db(10)=30: db(11)=45: db(12)=70
  ndata=11: np=500: lmda=632.99
  da(1, 0)=48.333:   da(1, 1)=21.055:  da(1, 2)=45
  da(2, 0)=51.658:   da(2, 1)=18.202:  da(2, 2)=47
  da(3, 0)=56.470:   da(3, 1)=15.337:  da(3, 2)=49
  da(4, 0)=65.146:   da(4, 1)=12.304:  da(4, 2)=51
  da(5, 0)=79.263:   da(5, 1)=10.238:  da(5, 2)=53
  da(6, 0)=99.569:   da(6, 1)=10.716:  da(6, 2)=55
  da(7, 0)=116.453:  da(7, 1)=13.864:  da(7, 2)=57
  da(8, 0)=128.599:  da(8, 1)=17.411:  da(8, 2)=59
  da(9, 0)=137.276:  da(9, 1)=20.076:  da(9, 2)=61
  da(10, 0)=145.740: da(10, 1)=21.970: da(10, 2)=65
  da(11, 0)=145.758: da(11, 1)=22.216: da(11, 2)=70
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b2]                                                   'Data2 parameterss
  nlay=2
  nsub=3.87: ksub=0.14: namb=1.00027
  db(1)=1.48: db(2)=0: db(3)=100: db(4)=1.58: db(5)=0: db(6)=100
  db(7)=45: db(8)=95: db(9)=10: db(10)=40: db(11)=45: db(12)=70
  ndata=6: np=500: lmda=632.99
  da(1, 0)=50.160:   da(1, 1)=34.744:  da(1, 2)=45
  da(2, 0)=53.459:   da(2, 1)=31.895:  da(2, 2)=50
  da(3, 0)=58.421:   da(3, 1)=28.871:  da(3, 2)=55
  da(4, 0)=65.798:   da(4, 1)=25.852:  da(4, 2)=60
  da(5, 0)=76.205:   da(5, 1)=23.761:  da(5, 2)=65
  da(6, 0)=89.334:   da(6, 1)=23.740:  da(6, 2)=70
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b3]
  if kk>0 and kk<4 then
    gosub [getnums]
    #1 "color darkpink": gosub [tplot]
  end if
goto [waitHere]

[b4]
  if kk>0 and kk<4 then
    gosub [getnums]
    kk=3: #1 "color green": gosub [tplot]
    goto [sim]
  end if
goto [waitHere]

[b5]
  if kk=3 then
    kk=4: r=105: c=90
    #1 "cls ; color black"
    #1 "place "; c; " "; r: #1 "\P  deg": c=c+120
    #1 "place "; c; " "; r: #1 "\A  deg": c=c+120
    #1 "place "; c; " "; r: #1 "\Theta deg": c=c+120
    #1 "place "; c; " "; r: #1 "\Dev  deg": c=90
    for i=1 to ndata
      r=r+15
      for j=0 to 3
        #1 "place "; c+j*120; " "; r: #1 "\"; da(i, j)
      next j
    next i
    r=r+15
    #1 "place 390 "; r: #1 "\rms dev"
    #1 "place 450 "; r: #1 "\"; q(1, sn)
goto [waitHere]

[grid]
  #1 "cls ; down ; color black"
  for i=0 to 5
    x2=90+i*100
    #1 "line "; x2; " 90 "; x2; " 390"
  next i
  for i=0 to 3
    y2=90+i*100
    #1 "line 90 "; y2; " 590 "; y2
  next i
  #1 "place 85 410":  #1 "\"; db(7)
  #1 "place 577 410": #1 "\"; db(8)
  #1 "place 67 396":  #1 "\"; db(9)
  #1 "place 67 96":   #1 "\"; db(10)
  #1 "place 320 410": #1 "\P  deg"
  #1 "place 60 238":  #1 "\A"
  #1 "place 54 258":  #1 "\deg"
  #1 "backcolor red"
  for i=1 to ndata
    pp=da(i, 0): aa=da(i,1): gosub [scale]
    #1 "place "; xs; " "; ys: #1 "circlefilled 5"
  next i
  #1 "backcolor white"
return

[getnums]
  for i=1 to sn: q(1, i)=0: next i
  #1.t1, "!contents?": input#1.t1, z$:  q(1, 1)=val(z$)
  #1.t2, "!contents?": input#1.t2, z$:  q(1, 2)=val(z$)
  #1.t3, "!contents?": input#1.t3, z$:  q(1, 3)=val(z$)
  if nlay=2 then
    #1.t4, "!contents?": input#1.t4, z$:  q(1, 4)=val(z$)
    #1.t5, "!contents?": input#1.t5, z$:  q(1, 5)=val(z$)
    #1.t6, "!contents?": input#1.t6, z$:  q(1, 6)=val(z$)
  end if
return

[setnums]
  #1.t1, str$(db(1)): #1.t2, str$(db(2)): #1.t3, str$(db(3))
  if nlay=2 then
    #1.t4, str$(db(4)): #1.t5, str$(db(5)): #1.t6, str$(db(6))
  else
    #1.t4, " ": #1.t5, " ": #1.t6, " "
  end if
return

[scale]
  xs=90+int((pp-db(7))*500/(db(8)-db(7)))
  ys=390-int((aa-db(9))*300/(db(10)-db(9)))
return

[tplot]
  pold=0: iv=1
  theta=db(11)
  gosub [pa]
  gosub [scale]: #1 "place "; xs; " "; ys
  for ii=1 to np
    theta=db(11)+ii*(db(12)-db(11))/np: gosub [pa]: gosub [scale]
    if abs(pp-pold)>120 then z$="place " else z$="goto "
    #1 z$; xs; " "; ys: pold=pp
  next ii
  #1 "color black ; backcolor yellow"
  for ii=1 to ndata
    theta=da(ii, 2): gosub [pa]: gosub [scale]
    #1 "place "; xs; " "; ys: #1 "circlefilled 5"
  next ii
  #1 "backcolor white"
return

[pa]                                                   ' Calculate P and A
  j(0)=nsub: k(0)=ksub
  for i=1 to nlay
    j(i)=q(iv, 3*i-2): k(i)=q(iv, 3*i-1)
  next i
  bb=namb*sin(theta*pi/180): bb=bb*bb
  b(0)=namb*cos(theta*pi/180): b(1)=b(0)/namb/namb
  for i=0 to nlay
    e2=j(i)*j(i): e3=k(i)*k(i): e4=e2-e3-bb
    e5=sqr(e4*e4+4*e2*e3): e6=1/(e2+e3)/(e2+e3)
    c(i, 0)=sqr((e4+e5)/2)
    d(i, 0)=0-j(i)*k(i)/c(i, 0)
    c(i, 1)=e6*(c(i, 0)*(e2-e3)-2*d(i, 0)*k(i)*j(i))
    d(i, 1)=e6*(d(i, 0)*(e2-e3)+2*c(i, 0)*k(i)*j(i))
  next i
  for i=0 to 1
    for il=1 to nlay
      f1=q(iv, 3*il)*2*pi/lmda: f2=f1*c(il, 0): f3=f1*d(il, 0)
      f4=sin(f2): f5=cos(f2): f6=exp(f3)/2
      f7=exp(0-f3)/2: e7=f5*(f7-f6): e8=f4*(f6+f7)
      if il>1 then
        for is=0 to 3
          m(is)=e(is): h(is)=f(is)
        next is
      end if
      e(0)=f5*(f6+f7): f(0)=f4*(f7-f6): f(2)=0-e7*c(il, i)+e8*d(il, i)
      e(2)=e7*d(il, i)+e8*c(il, i): e9=c(il, i)*c(il, i)+d(il, i)*d(il, i)
      f(1)=(e7*c(il, i)+e8*d(il, i))/e9: e(1)=(0-e8*c(il, i)+e7*d(il, i))/e9
      e(3)=e(0): f(3)=f(0)
      if il>1 then
        for im=0 to 1
          for in=0 to 2 step 2
            r(im+in)=e(in)*m(im)-f(in)*h(im)+e(in+1)*m(im+2)-f(in+1)*h(im+2)
            s(im+in)=f(in)*m(im)+e(in)*h(im)+e(in+1)*h(im+2)+f(in+1)*m(im+2)
          next in
        next im
        for is=0 to 3
          e(is)=r(is): f(is)=s(is)
        next is
      end if
    next il
    g1=(e(0)+f(1)*c(0, i)+e(1)*d(0, i))*b(i)
    g2=(f(0)+f(1)*d(0, i)-e(1)*c(0, i))*b(i)
    g3=0-f(2)+e(3)*c(0, i)-f(3)*d(0, i)
    g4=e(2)+e(3)*d(0, i)+f(3)*c(0, i)
    g5=(g1+g3)*(g1+g3)+(g2+g4)*(g2+g4)
    u(i)=(g1*g1+g2*g2-g3*g3-g4*g4)/g5: v(i)=2*(g2*g3-g1*g4)/g5
  next i
  g6=u(0)*u(0)+v(0)*v(0): g7=(u(0)*u(1)+v(0)*v(1))/g6
  g8=(u(0)*v(1)-u(1)*v(0))/g6
  aa=atn((g7*g7+g8*g8)^.5)*180/pi
  pp=atn(g8/g7)*180/pi: if g7<0 then pp=pp+180
  pp=pp/2-45
  if pp<db(7) then pp=pp+180
return

[quit]
  close #1
end