' simpani.bas - jack@ord.ca - revised 16 Feb 04 - Liberty Basic v3.01
' the algorithm optimizes the fit of P-A data to the growth curve for
'   a transparent anisotropic film
' NOTES: (1) THE ALGORITHM DOES ARITHMETIC, BUT MAKES NO CLAIM ABOUT
'            THE SIGNIFICANCE OF THE PARAMETER VALUES IT DETERMINES
'        (2) "Data2" REQUIRES THE FILE "Va3o.txt" IN THE SAME DIRECTORY
'        (3) "Print" LISTS ONLY EVERY THIRD POINT IN "Data2"

' 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, "Nx", [waitHere], [waitHere], 190, 3, 50, 20
  checkbox#1.c4, "Nz", [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, "Tplot", [b3], UL, 1, 116, 50, 20
  button#1, "Smplx", [b4], UL, 1, 140, 50, 20
  button#1, "Print", [b5], UL, 1, 164, 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 anisotropic films" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim ip(6): dim db(13)
  dim c(3, 2): dim d(3, 2)
  dim q(9, 6): dim da(70, 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
  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.372: db(0)=0
  db(1)=4.07996: db(2)=3.20502: db(3)=2.1692: db(4)=2.1692
  db(5)=0: db(6)=160
  db(7)=0: db(8)=80: db(9)=25: db(10)=75
  ndata=16: np=500: lmda=632.99
  da(1, 0)=24.708:   da(1, 1)=29.561                   ' P-A data
  da(2, 0)=19.052:   da(2, 1)=30.594
  da(3, 0)=14.164:   da(3, 1)=32.142
  da(4, 0)=10.341:   da(4, 1)=34.574
  da(5, 0)=8.148:    da(5, 1)=38.538
  da(6, 0)=8.608:    da(6, 1)=45.327
  da(7, 0)=14.639:   da(7, 1)=56.898
  da(8, 0)=37.582:   da(8, 1)=68.422
  da(9, 0)=65.106:   da(9, 1)=59.783
  da(10, 0)=72.847:  da(10, 1)=47.007
  da(11, 0)=73.697:  da(11, 1)=39.070
  da(12, 0)=71.540:  da(12, 1)=34.620
  da(13, 0)=67.584:  da(13, 1)=31.931
  da(14, 0)=62.531:  da(14, 1)=30.274
  da(15, 0)=56.571:  da(15, 1)=29.214
  da(16, 0)=49.834:  da(16, 1)=28.604
  kk=1: gosub [setnums]: gosub [grid]
goto [waitHere]

[b2]                                                   ' Data2 parameters
  theta=60: namb=1.372: db(0)=0
  db(1)=3.637: db(2)=3.334: db(3)=2.2902: db(4)=2.2902
  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]

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

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

[b5]
  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)
goto [waitHere]

[grid]
  #1 "cls": #1 "down": #1 "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 50 208":  #1 "\A"
  #1 "place 44 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
    n=q(iv, 1): k=q(iv, 2)
    e2=n*n: e3=k*k: e4=e2-e3-bb: e5=(e4*e4+4*e2*e3)^.5
    c(0, 0)=((e4+e5)/2)^.5: d(0, 0)=0-n*k/c(0, 0): e6=1/(e2+e3)/(e2+e3)
    c(0, 1)=e6*(c(0, 0)*(e2-e3)-2*d(0, 0)*n*k)
    d(0, 1)=e6*(d(0, 0)*(e2-e3)+2*c(0, 0)*n*k)
    nx=q(iv, 3): nz=q(iv, 4)
    c(1, 0)=(nx*nx-bb)^.5: d(1, 0)=c(1, 0)
    c(1, 1)=nx*(1-bb/nz/nz)^.5: d(1, 1)=c(1, 1)/nx/nx
return

[pa]                                                   ' Calculate P and A
  for i=0 to 1
    f1=l*2*pi/lmda*c(1, i): f2=sin(f1): f3=cos(f1)
    e(0)=f3: e(1)=0-f2/d(1, i): e(2)=f2*d(1, i): e(3)=e(0)
    g1=(e(0)+e(1)*d(0, i))*b(i)
    g2=0-e(1)*c(0, i)*b(i)
    g3=e(3)*c(0, i)
    g4=e(2)+e(3)*d(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 "Va3o.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