' simplex.bas - jack@ord.ca - revised 11 Apr 2011 - Liberty Basic v3.01
' the algorithm optimizes the fit of P-A data to the growth curve for
'   a single layer
' NOTES: (1) THE ALGORITHM DOES ARITHMETIC, BUT MAKES NO CLAIM ABOUT
'            THE SIGNIFICANCE OF THE PARAMETER VALUES IT DETERMINES
'        (2) "Data3" REQUIRES THE FILE "Va3r.txt" IN THE SAME DIRECTORY
'        (3) "Print" LISTS ONLY EVERY THIRD POINT IN "Data3"

' Initialize Window
  nomainwin
  textbox#1.t1, 51, 3, 100, 20:   textbox#1.t2, 51, 27, 100, 20
  textbox#1.t3, 240, 3, 100, 20:  textbox#1.t4, 240, 27, 100, 20
  checkbox#1.c1, "Ns", [waitHere], [waitHere], 1, 3, 50, 20
  checkbox#1.c2, "Ks", [waitHere], [waitHere], 1, 27, 50, 20
  checkbox#1.c3, "Nf", [waitHere], [waitHere], 190, 3, 50, 20
  checkbox#1.c4, "Kf", [waitHere], [waitHere], 190, 27, 50, 20
  button#1, "Data1", [b1], UL, 1, 68, 50, 20
  button#1, "Data2", [b2], UL, 1, 92, 50, 20
  button#1, "Data3", [b3], UL, 1, 116, 50, 20
  button#1, "Tplot", [b4], UL, 1, 140, 50, 20
  button#1, "Smplx", [b5], UL, 1, 164, 50, 20
  button#1, "Print", [b6], UL, 1, 188, 50, 20
  graphicbox#1.g1, 380, 1, 200, 50
  WindowWidth=610     ' pixel scale 0-600
  WindowHeight=429    ' pixel scale 0-400
  UpperLeftX=10: UpperLeftY=10
  open "Simplex algorithm for film-growth data" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim ip(6): dim db(13)
  dim b(2): dim c(3, 2): dim d(3, 2)
  dim q(9, 6): dim da(55, 4)
  dim u(2): dim v(2): dim e(4): dim f(4)
  pi=4*atn(1): kk=0: w=0: sn=5: f$="##.######"
  #1 "down ; place 80 60 ; box 580 360"

[waitHere]
  wait

[sim]
  #1.g1 "cls ; down ; color black"
  for i=1 to 3
    x1=50*i: #1.g1 "line "; x1; " 0 "; x1; " 50"
  next i
  #1.g1 "color red"
  i=0: jlof=0: test=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
  #1.c4, "value?": input#1.c4, z$: if z$="set" then i=i+1: ip(i)=4
  nvr=i: nvx=i+1: ip(nvx)=sn
  if nvr=0 then
    df=0: iv=1: gosub [lfit]
  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
      q(i, ip(i))=q(i, ip(i))+dni
    next i
    for iv=1 to nvx
      gosub [lfit]
    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=49 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 [lfit]
      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 [lfit]
            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 380": #1 "\rms dev   "; using(f$, q(1, sn))
  et=time$("ms" )
  #1 "place 400 380": #1 "\exec time "; et-tt; " msec"
  gosub [curve]
  #1.t1, using(f$, q(1, 1))
  #1.t2, using(f$, q(1, 2))
  #1.t3, using(f$, q(1, 3))
  #1.t4, using(f$, q(1, 4))
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 [lfit]
  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 [lfit]
    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

[lfit]
  dl=1: fs=0
  gosub [pq]
  for j=1 to ndata
    la=da(j, 2)
    if (la=0) and (j>1) then la=da(j-1, 2)
    l=la: gosub [dev]: fa=ds
    lb=la+dl: l=lb: gosub [dev]: fb=ds
    if fb>fa then lx=la: la=lb: lb=lx: fx=fa: fa=fb: fb=fx: dl=0-dl
    lx=lb+dl: l=lx: gosub [dev]: fx=ds
    while fx<fb
      lb=lx: fb=fx
      lx=lb+dl: l=lx: gosub [dev]: fx=ds
    wend
    cf=3: ddl=dl/cf: ddlf=0.01/cf
    while abs(ddl)>ddlf
      lx=lb+ddl: l=lx: gosub [dev]: fx=ds
      if fx<fb then
        lb=lx: fb=fx
      else
        ddl=0-ddl/cf
      end if
    wend
    da(j, 2)=lb: da(j, 3)=fb^.5: fs=fs+fb
  next j
  q(iv, sn)=(fs/ndata)^.5
return

[b1]                                                   ' Data1 parameters
  theta=60: namb=1.00027: nlay=1: db(0)=0
  db(1)=3.875: db(2)=.023: db(3)=1.75: db(4)=0: db(5)=0: db(6)=180
  db(7)=0-45: db(8)=135: db(9)=0: db(10)=90
  ndata=18: np=500: lmda=632.99
  da(1, 0)=33.418:    da(1, 1)=23.833                   ' P-A data
  da(2, 0)=28.385:    da(2, 1)=24.082
  da(3, 0)=22.334:    da(3, 1)=24.848
  da(4, 0)=18.46:     da(4, 1)=25.474
  da(5, 0)=14.078:    da(5, 1)=26.586
  da(6, 0)=8.916:     da(6, 1)=28.476
  da(7, 0)=1.057:     da(7, 1)=34.151
  da(8, 0)=0-8.117:   da(8, 1)=45.118
  da(9, 0)=0-11.112:  da(9, 1)=49.126
  da(10, 0)=0-15.396: da(10, 1)=53.622
  da(11, 0)=0-31.28:  da(11, 1)=63.657
  da(12, 0)=122.834:  da(12, 1)=64.013
  da(13, 0)=108.742:  da(13, 1)=56.62
  da(14, 0)=105.317:  da(14, 1)=53.373
  da(15, 0)=93.518:   da(15, 1)=39.348
  da(16, 0)=85.66:    da(16, 1)=31.144
  da(17, 0)=79.404:   da(17, 1)=27.602
  da(18, 0)=70.383:   da(18, 1)=24.847
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b2]                                                   'Data2 parameterss
  theta=60: namb=1.3341: nlay=1: db(0)=0
  db(1)=2.39: db(2)=4.236: db(3)=1.5: db(4)=.3: db(5)=0: db(6)=300
  db(7)=0-60: db(8)=30: db(9)=0: db(10)=60
  ndata=18: np=500: lmda=632.99
  da(1, 0)=17.779:    da(1, 1)=32.529                 ' P-A data
  da(2, 0)=11.909:    da(2, 1)=31.462
  da(3, 0)=7.456:     da(3, 1)=32.056
  da(4, 0)=3.791:     da(4, 1)=34.147
  da(5, 0)=.006:      da(5, 1)=37.689
  da(6, 0)=0-5.038:   da(6, 1)=42.459
  da(7, 0)=0-12.612:  da(7, 1)=47.444
  da(8, 0)=0-23.629:  da(8, 1)=50.099
  da(9, 0)=0-36.234:  da(9, 1)=47.329
  da(10, 0)=0-46.105: da(10, 1)=39.356
  da(11, 0)=0-51.844: da(11, 1)=29.789
  da(12, 0)=0-53.83:  da(12, 1)=21.196
  da(13, 0)=0-51.869: da(13, 1)=14.641
  da(14, 0)=0-44.808: da(14, 1)=10.464
  da(15, 0)=0-33.738: da(15, 1)=8.972
  da(16, 0)=0-24.041: da(16, 1)=9.634
  da(17, 0)=0-18.659: da(17, 1)=11.271
  da(18, 0)=0-16.169: da(18, 1)=13.047
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b3]                                                   ' Data3 parameters
  theta=60: namb=1.372: nlay=1: db(0)=0
  db(1)=3.637: db(2)=3.334: db(3)=2.0: db(4)=0
  db(5)=0: db(6)=120
  db(7)=0: db(8)=80: db(9)=25: db(10)=75
  np=500: lmda=632.99
  gosub [getdata]
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b4]
  if kk=1 or kk=2 then
    gosub [getnums]: #1 "color darkpink": gosub [curve]
  end if
goto [waitHere]

[b5]
  if kk=1 or kk=2 then
    kk=2: gosub [getnums]
    #1 "color green": gosub [curve]
    goto [sim]
  end if
goto [waitHere]

[b6]
  if kk=2 then
    kk=3: r=60: 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 "\Thick  nm": c=c+120
    #1 "place "; c; " "; r: #1 "\Dev  deg": c=90: r=65
    if ndata>20 then ii=3 else ii=1
    for i=1 to ndata step ii
      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)
  end if
goto [waitHere]

[grid]
  #1 "cls ; down ; color black"
  for i=0 to 5
    x2=80+i*100
    #1 "line "; x2; " 60 "; x2; " 360"
  next i
  for i=0 to 3
    y2=60+i*100
    #1 "line 80 "; y2; " 580 "; y2
  next i
  #1 "place 70 380":  #1 "\"; db(7)
  #1 "place 570 380": #1 "\"; db(8)
  #1 "place 60 366":  #1 "\"; db(9)
  #1 "place 60 66":   #1 "\"; db(10)
  #1 "place 310 380": #1 "\P  deg"
  #1 "place 54 208":  #1 "\A"
  #1 "place 46 228":  #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]
  #1.t1,  "!contents?": input#1.t1, z$:  db(1)=val(z$)
  #1.t2,  "!contents?": input#1.t2, z$:  db(2)=val(z$)
  #1.t3,  "!contents?": input#1.t3, z$:  db(3)=val(z$)
  #1.t4,  "!contents?": input#1.t4, z$:  db(4)=val(z$)
  for i=1 to 4
    q(1, i)=db(i)
  next i
  for i=1 to ndata
    da(i, 2)=0
  next i
return

[setnums]
  b1=theta*pi/180: bb=namb*sin(b1): bb=bb*bb
  b(0)=namb*cos(b1): b(1)=b(0)/namb/namb
  if db(0)>0 then
    b2=(2*db(1)+90)*pi/180: b3=tan(db(2)*pi/180)
    c1=bb/b(0)/(b3*b3+2*b3*cos(b2)+1): c2=c1*(1-b3*b3)
    c3=2*b3*c1*sin(b2): c4=c2*c2-c3*c3+bb
    c5=2*c2*c3: c6=(c4*c4+c5*c5)^.5
    db(1)=((c4+c6)/2)^.5: db(2)=c5/2/db(1)
  end if
  #1.t1, str$(db(1))
  #1.t2, str$(db(2))
  #1.t3, str$(db(3))
  #1.t4, str$(db(4))
return

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

[curve]
  pold=0: iv=1: gosub [pq]
  l=db(5): gosub [pa]: gosub [scale]: #1 "place "; xs; " "; ys
  for ii=1 to np
    l=db(5)+ii*(db(6)-db(5))/np: gosub [pa]: gosub [scale]
    if abs(pp-pold)>90 then z$="place " else z$="goto "
    #1 z$; xs; " "; ys: pold=pp
  next ii
return

[pq]                                                   ' Calculate p and q
  for i=0 to nlay
    n=q(iv, 2*i+1): k=q(iv, 2*i+2)
    e2=n*n: e3=k*k: e4=e2-e3-bb: e5=(e4*e4+4*e2*e3)^.5
    c(i, 0)=((e4+e5)/2)^.5: d(i, 0)=0-n*k/c(i, 0): e6=1/(e2+e3)/(e2+e3)
    c(i, 1)=e6*(c(i, 0)*(e2-e3)-2*d(i, 0)*n*k)
    d(i, 1)=e6*(d(i, 0)*(e2-e3)+2*c(i, 0)*n*k)
  next i
return

[pa]                                                   ' Calculate P and A
  il=1
  f1=l*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)
  for i=0 to 1
    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)
    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

[getdata]
  open "Va3r.txt" for random as #2 len=12             ' get data from file
  ndata=lof(#2)/12
  field#2, 6 as a, 6 as b
  for i=1 to ndata
    get#2, i: da(i, 0)=a: da(i, 1)=b                  ' store in da(i, j)
  next i
  close#2
return

[quit]
  close #1
end