fini
/clear,start
*get,zeit,active,0,time,cpu
/DSCALE,ALL,1.0 
/config,nres,10000
save

!******************************
!Parameter input
!******************************
!geometry data
!******************************

!Number of finger electrodes
fingernum = 16

!Wavelength
wl = 0.56e-3

!width of electrodes (y direction)
ebreite = 1e-4

!substrate thickness
subdick = 3e-3

!Width of IDTs (z direction - for impedance determination)
b_IDT = 26.38e-3

!outer region for damping the SAW
prdaempf = 3*wl

!undamped traveling distance above actuator
laufstr = 4*wl

!undamped traveling distance under actuator
laufstr1 = wl/10

!thickness of top electrode
ded = 8e-6

!thickness of ground electrode
ged = 12e-6

!thickness of isolation layer
id = 25e-6

!thickness including substrate, isolation layer, ground electrode and PZT layer
kerdick = subdick+ged+id+1.51e-4

!overhang of PZT layer (y-direction)
pztueber = 0.16e-3

!overhang of isolation layer (y-direction)
isoueber = 1e-3

!******************************
!material data
!******************************

!Young's modulus substrate (DIN 1.4016)
emod = 2.1e11
!Poisson's ratio substrate
pz = 0.27
!Density substrate
dichte = 7850
!damping coefficient of outer region
daempf = 5e-9

!elasticity of PZT-8 (here given for Z-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
sz11 = 1.15e-11
sz12 = -3.70e-12
sz13 = -4.80e-12
sz14 = 0
sz15 = 0
sz16 = 0
sz21 = -3.70e-12
sz22 = 1.15e-11
sz23 = -4.80e-12
sz24 = 0
sz25 = 0
sz26 = 0
sz31 = -4.80e-12
sz32 = -4.80e-12
sz33 = 1.35e-11
sz34 = 0
sz35 = 0
sz36 = 0
sz41 = 0
sz42 = 0
sz43 = 0
sz44 = 3.19e-11
sz45 = 0
sz46 = 0
sz51 = 0
sz52 = 0
sz53 = 0
sz54 = 0
sz55 = 3.19e-11
sz56 = 0
sz61 = 0
sz62 = 0
sz63 = 0
sz64 = 0
sz65 = 0
sz66 = 3.04e-11
!piez. coupling-coeff. of PZT-8 (here given for Z-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
dz11 = 0
dz12 = 0
dz13 = -0.97e-10/225*140
dz21 = 0
dz22 = 0
dz23 = -0.97e-10/225*140
dz31 = 0
dz32 = 0
dz33 = 2.25e-10/225*140
dz41 = 0
dz42 = 3.30e-10/225*140
dz43 = 0
dz51 = 3.30e-10/225*140
dz52 = 0
dz53 = 0
dz61 = 0
dz62 = 0
dz63 = 0
!rel. dielektr. const. of PZT-8 (here given for Z-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
erz11 = 1290
erz12 = 0
erz13 = 0
erz21 = 0
erz22 = 1290
erz23 = 0
erz31 = 0
erz32 = 0
erz33 = 1000

!elasticity of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
hsx11 = sz33
hsx12 = sz31
hsx13 = sz32
hsx14 = sz36
hsx15 = sz34
hsx16 = sz35
hsx21 = sz13
hsx22 = sz11
hsx23 = sz12
hsx24 = sz16
hsx25 = sz14
hsx26 = sz15
hsx31 = sz23
hsx32 = sz21
hsx33 = sz22
hsx34 = sz26
hsx35 = sz24
hsx36 = sz25
hsx41 = sz63
hsx42 = sz61
hsx43 = sz62
hsx44 = sz66
hsx45 = sz64
hsx46 = sz65
hsx51 = sz43
hsx52 = sz41
hsx53 = sz42
hsx54 = sz46
hsx55 = sz44
hsx56 = sz45
hsx61 = sz53
hsx62 = sz51
hsx63 = sz52
hsx64 = sz56
hsx65 = sz54
hsx66 = sz55
!piez. coupling-coeff. of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
hdx11 = dz33
hdx12 = dz31
hdx13 = dz32
hdx21 = dz13
hdx22 = dz11
hdx23 = dz12
hdx31 = dz23
hdx32 = dz21
hdx33 = dz22
hdx41 = dz63
hdx42 = dz61
hdx43 = dz62
hdx51 = dz43
hdx52 = dz41
hdx53 = dz42
hdx61 = dz53
hdx62 = dz51
hdx63 = dz52
!rel. dielektr. const. of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around x,5..around y,6..around z)
herx11 = erz33
herx12 = erz31
herx13 = erz32
herx21 = erz13
herx22 = erz11
herx23 = erz12
herx31 = erz23
herx32 = erz21
herx33 = erz22

!elasticity of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around z,5..around x,6..around y)
sx11 = hsx11
sx12 = hsx12
sx13 = hsx13
sx14 = hsx16
sx15 = hsx14
sx16 = hsx15
sx21 = hsx21
sx22 = hsx22
sx23 = hsx23
sx24 = hsx26
sx25 = hsx24
sx26 = hsx25
sx31 = hsx31
sx32 = hsx32
sx33 = hsx33
sx34 = hsx36
sx35 = hsx34
sx36 = hsx35
sx41 = hsx61
sx42 = hsx62
sx43 = hsx63
sx44 = hsx66
sx45 = hsx64
sx46 = hsx65
sx51 = hsx41
sx52 = hsx42
sx53 = hsx43
sx54 = hsx46
sx55 = hsx44
sx56 = hsx45
sx61 = hsx51
sx62 = hsx52
sx63 = hsx53
sx64 = hsx56
sx65 = hsx54
sx66 = hsx55
!piez. coupling-coeff. of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around z,5..around x,6..around y)
dx11 = hdx11
dx12 = hdx12
dx13 = hdx13
dx21 = hdx21
dx22 = hdx22
dx23 = hdx23
dx31 = hdx31
dx32 = hdx32
dx33 = hdx33
dx41 = hdx61
dx42 = hdx62
dx43 = hdx63
dx51 = hdx41
dx52 = hdx42
dx53 = hdx43
dx61 = hdx51
dx62 = hdx52
dx63 = hdx53
!rel. dielektr. const. of PZT-8 (here given for X-polarization; format: 1..x,2..y,3..z,4..around z,5..around x,6..around y)
erx11 = herx11
erx12 = herx12
erx13 = herx13
erx21 = herx21
erx22 = herx22
erx23 = herx23
erx31 = herx31
erx32 = herx32
erx33 = herx33

!Density of PZT-8
dichtepzt = 7600

!Young's modulus of elektrode (gold)
emodele = 7e10
!Poisson's ratio of electrode
pzele = 0.42
!Density of electrode
dichteele = 19320

!Young's modulus of isolation (Macor)
emodiso = 6.69e10
!Poisson's ratio of isolation
pziso = 0.29
!Density of isolation
dichteiso = 2520
!rel. dielektr. const. of isolation
erz11iso = 5


!******************************
!stimulation and simulation data
!******************************
!Voltage amplitude (0-peak)
u0p = 1

!Start frequency
fmin = 3.5e6

!End frequency
fmax = 6.0e6

!Number of series points
zyxend=50

!Number of time periods per series point
nper = 40

!Number of periods, considered for results
AnzMW = 10



!******************************
!Creating output file
!******************************
*get,ctime1,active,0,time,cpu	!generating unique ID
*get,ctime,active,0,time,wall
ctime=ctime+ctime1
/OUTPUT,_output_%ctime%,out
*cfopen,Deflections_%fingernum%_%wl%_%ctime%,txt
*cfwrite,,Frequency,Num_fingers,average_ampl,standard_derivation_Ampl,Re_Z,Im_Z
*cfclos
save


*do,zyx,0,zyxend	!Calculate each series point
 parsav,all
 resume
 parres
 
 /prep7
 
 
 *if,zyx,eq,0,then	!Defining arrays for calculation of averaged deflection
  *DIM,disknot,array,nint(laufstr/wl-1),AnzMW
  *DIM,knotnumarr,array,nint(laufstr/wl-1)
  *DIM,flagzaehl,array,nint(laufstr/wl-1),2
  *DIM,maxux,array,nint(laufstr/wl-1)
 *endif
 
 
 !******************************
 !resulting parameters
 !******************************
 
 !Pi
 pi=acos(-1)
 
 !current frequency
 f = fmin+(fmax-fmin)*zyx/zyxend
 
 !mesh size
 vdichte = wl/10/1.25**4
 
 !simulation time
 tsim = nper/f
 
 !simulation steps
 schritte = 40*tsim*f
 
 !angular velocity
 omega = 2*pi*f
 
 !******************************
 !geometry definition
 !******************************
 !key points
 !******************************
 
 k,1,0,0,0
 k,2,subdick,0,0
 
 k,3,0,prdaempf,0
 k,4,subdick,prdaempf,0
 
 k,5,0,prdaempf+laufstr1,0
 k,6,subdick,prdaempf+laufstr1,0
 k,7,subdick+id,prdaempf+laufstr1,0
 k,8,subdick+id+ged,prdaempf+laufstr1,0
 
 k,9,0,prdaempf+laufstr1+isoueber,0
 k,10,subdick,prdaempf+laufstr1+isoueber,0
 k,11,subdick+id,prdaempf+laufstr1+isoueber,0
 k,12,subdick+id+ged,prdaempf+laufstr1+isoueber,0
 k,13,kerdick,prdaempf+laufstr1+isoueber,0
 
 *do,i,1,fingernum
  k,12*i+2,0,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+3,subdick,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+4,subdick+id,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+5,subdick+id+ged,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+6,kerdick,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+7,kerdick+ded,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2,0
  k,12*i+8,0,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
  k,12*i+9,subdick,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
  k,12*i+10,subdick+id,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
  k,12*i+11,subdick+id+ged,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
  k,12*i+12,kerdick,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
  k,12*i+13,kerdick+ded,prdaempf+laufstr1+isoueber+pztueber+(i-1)*wl/2+ebreite,0
 *enddo
 
 yhilf = prdaempf+laufstr1+isoueber+pztueber+(fingernum-1)*wl/2+ebreite+pztueber
 
 k,12*(fingernum)+14,0,yhilf,0
 k,12*(fingernum)+15,subdick,yhilf,0
 k,12*(fingernum)+16,subdick+id,yhilf,0
 k,12*(fingernum)+17,subdick+id+ged,yhilf,0
 k,12*(fingernum)+18,kerdick,yhilf,0
 k,12*(fingernum)+19,kerdick+ded,yhilf,0
 
 k,12*(fingernum)+20,0,yhilf+isoueber,0
 k,12*(fingernum)+21,subdick,yhilf+isoueber,0
 k,12*(fingernum)+22,subdick+id,yhilf+isoueber,0
 k,12*(fingernum)+23,subdick+id+ged,yhilf+isoueber,0
 
 k,12*(fingernum)+24,0,yhilf+isoueber+laufstr,0
 k,12*(fingernum)+25,subdick,yhilf+isoueber+laufstr,0
 
 k,12*(fingernum)+26,0,yhilf+isoueber+laufstr+prdaempf,0
 k,12*(fingernum)+27,subdick,yhilf+isoueber+laufstr+prdaempf,0
 
 !******************************
 !Areas
 !******************************
 a,1,2,4,3
 a,3,4,6,5
 a,5,6,10,9
 a,6,7,11,10
 a,7,8,12,11
 a,9,10,15,14
 a,10,11,16,15
 a,11,12,17,16
 a,12,13,18,17
 *do,i,1,fingernum
  a,12*i+2,12*i+3,12*i+9,12*i+8
  a,12*i+3,12*i+4,12*i+10,12*i+9
  a,12*i+4,12*i+5,12*i+11,12*i+10
  a,12*i+5,12*i+6,12*i+12,12*i+11
  a,12*i+6,12*i+7,12*i+13,12*i+12
  a,12*i+8,12*i+9,12*i+15,12*i+14
  a,12*i+9,12*i+10,12*i+16,12*i+15
  a,12*i+10,12*i+11,12*i+17,12*i+16
  a,12*i+11,12*i+12,12*i+18,12*i+17
 *enddo
 
 
 a,12*fingernum+14,12*fingernum+15,12*fingernum+21,12*fingernum+20
 a,12*fingernum+15,12*fingernum+16,12*fingernum+22,12*fingernum+21
 a,12*fingernum+16,12*fingernum+17,12*fingernum+23,12*fingernum+22
 a,12*fingernum+20,12*fingernum+21,12*fingernum+25,12*fingernum+24
 a,12*fingernum+24,12*fingernum+25,12*fingernum+27,12*fingernum+26
 
 
 !******************************
 !Determining mesh size
 !******************************
 locy = prdaempf/2
 lsel,s,loc,y,locy*(1-1e-5),locy*(1+1e-5)
 lsel,r,loc,x,0
 cm,daempflinien1,line	!defining components for separate mesh sizes
 lesize,all,vdichte*5,,,10
 locy = yhilf+isoueber+laufstr+prdaempf/2
 lsel,s,loc,y,locy*(1-1e-5),locy*(1+1e-5)
 lsel,r,loc,x,0
 cm,daempflinien2,line
 lesize,all,vdichte*5,,,1/10
 locy = prdaempf/2
 lsel,s,loc,y,locy*(1-1e-5),locy*(1+1e-5)
 lsel,r,loc,x,subdick
 cm,daempflinien3,line
 lesize,all,vdichte*5,,,1/10
 locy = yhilf+isoueber+laufstr+prdaempf/2
 lsel,s,loc,y,locy*(1-1e-5),locy*(1+1e-5)
 lsel,r,loc,x,subdick
 cm,daempflinien4,line
 lesize,all,vdichte*5,,,10
 
 lsel,s,loc,x,subdick/2
 lsel,u,,,1
 cm,tiefenlinien1,line
 lesize,all,vdichte*5,,,10
 lsel,s,,,1
 cm,tiefenlinien2,line
 lesize,all,vdichte*5,,,1/10
 
 allsel
 lsel,u,,,daempflinien1
 lsel,u,,,daempflinien2
 lsel,u,,,daempflinien3
 lsel,u,,,daempflinien4
 lsel,u,,,tiefenlinien1
 lsel,u,,,tiefenlinien2
 lesize,all,vdichte
 allsel
 
 !******************************
 !Defining components for different parts
 !******************************
 asel,s,loc,y,prdaempf/2
 asel,a,loc,y,yhilf+isoueber+laufstr+prdaempf/2
 cm,Daempfung,area
 asel,s,loc,x,subdick/2
 asel,u,,,Daempfung
 cm,substrate,area
 asel,s,loc,x,subdick+id/2
 cm,Isolation,area
 asel,s,loc,x,subdick+id+ged/2
 cm,Leitschicht,area
 asel,s,loc,x,subdick+id+ged+(kerdick-subdick-id-ged)/2
 cm,Piezo,area
 asel,s,loc,x,14,kerdick+ded/2
 cm,Elektr,area
 
 allsel
 
 
 !******************************
 !Allocate and mesh substrate
 !******************************
 et,1,223
 keyopt,1,1,1001
 keyopt,1,3,2
 asel,s,,,substrate
 
 aktemod = emod
 aktpz = pz
 aktdichte = dichte
 
 TB,ANEL,3,1,21,1
 TBTEMP,0
 TBDATA,,1/aktemod,-1*aktpz/aktemod,-1*aktpz/aktemod,0,0,0  
 TBDATA,,1/aktemod,-1*aktpz/aktemod,0,0,0,1/aktemod   
 TBDATA,,0,0,0,2*(1+aktpz)/aktemod,0,0  
 TBDATA,,2*(1+aktpz)/aktemod,0,2*(1+aktpz)/aktemod,,, 
 TBDE,PIEZ,1 
 TB,PIEZ,3,,,1   
 TBMODIF,1,1,0
 TBMODIF,1,2,0
 TBMODIF,1,3,0
 TBMODIF,2,1,0
 TBMODIF,2,2,0
 TBMODIF,2,3,0
 TBMODIF,3,1,0
 TBMODIF,3,2,0
 TBMODIF,3,3,0
 TBMODIF,4,1,0
 TBMODIF,4,2,0
 TBMODIF,4,3,0
 TBMODIF,5,1,0
 TBMODIF,5,2,0
 TBMODIF,5,3,0
 TBMODIF,6,1,0
 TBMODIF,6,2,0
 TBMODIF,6,3,0
 TBDE,DPER,1 
 TB,DPER,3,,,1   
 TBMODIF,1,1,1
 TBMODIF,1,2,1
 TBMODIF,1,3,1
 TBMODIF,1,4,0   
 TBMODIF,1,5,0   
 TBMODIF,1,6,0    
 MP,DENS,3,aktdichte
 
 type,1
 mat,3
 real,1
 amesh,all
 
 !******************************
 !Allocate and mesh outer regions of substrate
 !******************************
 et,4,223
 keyopt,4,1,1001
 keyopt,4,3,2
 asel,s,,,Daempfung
 
 aktemod = emod
 aktpz = pz
 aktdichte = dichte
 
 TB,ANEL,4,1,21,1
 TBTEMP,0
 TBDATA,,1/aktemod,-1*aktpz/aktemod,-1*aktpz/aktemod,0,0,0  
 TBDATA,,1/aktemod,-1*aktpz/aktemod,0,0,0,1/aktemod   
 TBDATA,,0,0,0,2*(1+aktpz)/aktemod,0,0  
 TBDATA,,2*(1+aktpz)/aktemod,0,2*(1+aktpz)/aktemod,,, 
 TBDE,PIEZ,1 
 TB,PIEZ,4,,,1   
 TBMODIF,1,1,0
 TBMODIF,1,2,0
 TBMODIF,1,3,0
 TBMODIF,2,1,0
 TBMODIF,2,2,0
 TBMODIF,2,3,0
 TBMODIF,3,1,0
 TBMODIF,3,2,0
 TBMODIF,3,3,0
 TBMODIF,4,1,0
 TBMODIF,4,2,0
 TBMODIF,4,3,0
 TBMODIF,5,1,0
 TBMODIF,5,2,0
 TBMODIF,5,3,0
 TBMODIF,6,1,0
 TBMODIF,6,2,0
 TBMODIF,6,3,0
 TBDE,DPER,1 
 TB,DPER,4,,,1   
 TBMODIF,1,1,1
 TBMODIF,1,2,1
 TBMODIF,1,3,1
 TBMODIF,1,4,0   
 TBMODIF,1,5,0   
 TBMODIF,1,6,0    
 MP,DENS,4,aktdichte
 mp,damp,4,daempf
 
 type,4
 mat,4
 real,1
 amesh,all
 
 !******************************
 !Allocate and mesh electrodes
 !******************************
 et,5,223
 keyopt,5,1,1001
 keyopt,5,3,2
 asel,s,,,Elektr
 asel,a,,,Leitschicht
 
 aktemod = emodele
 aktpz = pzele
 aktdichte = dichteele
 
 TB,ANEL,5,1,21,1
 TBTEMP,0
 TBDATA,,1/aktemod,-1*aktpz/aktemod,-1*aktpz/aktemod,0,0,0  
 TBDATA,,1/aktemod,-1*aktpz/aktemod,0,0,0,1/aktemod   
 TBDATA,,0,0,0,2*(1+aktpz)/aktemod,0,0  
 TBDATA,,2*(1+aktpz)/aktemod,0,2*(1+aktpz)/aktemod,,, 
 TBDE,PIEZ,1 
 TB,PIEZ,5,,,1   
 TBMODIF,1,1,0
 TBMODIF,1,2,0
 TBMODIF,1,3,0
 TBMODIF,2,1,0
 TBMODIF,2,2,0
 TBMODIF,2,3,0
 TBMODIF,3,1,0
 TBMODIF,3,2,0
 TBMODIF,3,3,0
 TBMODIF,4,1,0
 TBMODIF,4,2,0
 TBMODIF,4,3,0
 TBMODIF,5,1,0
 TBMODIF,5,2,0
 TBMODIF,5,3,0
 TBMODIF,6,1,0
 TBMODIF,6,2,0
 TBMODIF,6,3,0
 TBDE,DPER,1 
 TB,DPER,5,,,1   
 TBMODIF,1,1,1
 TBMODIF,1,2,1
 TBMODIF,1,3,1
 TBMODIF,1,4,0   
 TBMODIF,1,5,0   
 TBMODIF,1,6,0    
 MP,DENS,5,aktdichte
 
 !******************************
 !Allocate and mesh isolation
 !******************************
 type,5
 mat,5
 real,4
 amesh,all
 
 et,6,223
 keyopt,6,1,1001
 keyopt,6,3,2
 asel,s,,,Isolation
 
 aktemod = emodiso
 aktpz = pziso
 aktdichte = dichteiso
 
 TB,ANEL,6,1,21,1
 TBTEMP,0
 TBDATA,,1/aktemod,-1*aktpz/aktemod,-1*aktpz/aktemod,0,0,0  
 TBDATA,,1/aktemod,-1*aktpz/aktemod,0,0,0,1/aktemod   
 TBDATA,,0,0,0,2*(1+aktpz)/aktemod,0,0  
 TBDATA,,2*(1+aktpz)/aktemod,0,2*(1+aktpz)/aktemod,,, 
 TBDE,PIEZ,1 
 TB,PIEZ,6,,,1   
 TBMODIF,1,1,0!1e-30
 TBMODIF,1,2,0
 TBMODIF,1,3,0
 TBMODIF,2,1,0!-1e-30
 TBMODIF,2,2,0
 TBMODIF,2,3,0
 TBMODIF,3,1,0!-1e-30
 TBMODIF,3,2,0
 TBMODIF,3,3,0
 TBMODIF,4,1,0
 TBMODIF,4,2,0
 TBMODIF,4,3,0
 TBMODIF,5,1,0
 TBMODIF,5,2,0
 TBMODIF,5,3,0!1e-30
 TBMODIF,6,1,0
 TBMODIF,6,2,0!1e-30
 TBMODIF,6,3,0
 TBDE,DPER,1 
 TB,DPER,6,,,1   
 TBMODIF,1,1,erz11iso
 TBMODIF,1,2,erz11iso
 TBMODIF,1,3,erz11iso
 TBMODIF,1,4,0   
 TBMODIF,1,5,0   
 TBMODIF,1,6,0    
 MP,DENS,6,aktdichte
 
 type,6
 mat,6
 real,6
 amesh,all
 
 !******************************
 !Allocate and mesh PZT (polarized in X-direction)
 !******************************
 asel,s,,,piezo
 et,2,223
 keyopt,2,1,1001
 keyopt,2,3,2
 TB,ANEL,2,1,21,1
 TBTEMP,0
 TBDATA,,sx11,sx12,sx13,sx14,sx15,sx16  
 TBDATA,,sx22,sx23,sx24,sx25,sx26,sx33
 TBDATA,,sx34,sx35,sx36,sx44,sx45,sx46  
 TBDATA,,sx55,sx56,sx66,,, 
 TBDE,PIEZ,1 
 TB,PIEZ,2,,,1   
 TBMODIF,1,1,dx11
 TBMODIF,1,2,dx12
 TBMODIF,1,3,dx13
 TBMODIF,2,1,dx21
 TBMODIF,2,2,dx22
 TBMODIF,2,3,dx23
 TBMODIF,3,1,dx31
 TBMODIF,3,2,dx32
 TBMODIF,3,3,dx33
 TBMODIF,4,1,dx41
 TBMODIF,4,2,dx42
 TBMODIF,4,3,dx43
 TBMODIF,5,1,dx51
 TBMODIF,5,2,dx52
 TBMODIF,5,3,dx53
 TBMODIF,6,1,dx61
 TBMODIF,6,2,dx62
 TBMODIF,6,3,dx63
 TBDE,DPER,1 
 TB,DPER,2,,,1   
 TBMODIF,1,1,erx11	!relative permittivity required
 TBMODIF,1,2,erx22	!relative permittivity required
 TBMODIF,1,3,erx33	!relative permittivity required
 TBMODIF,1,4,erx12	!relative permittivity required
 TBMODIF,1,5,erx23	!relative permittivity required
 TBMODIF,1,6,erx13	!relative permittivity required
 MP,DENS,2,dichtepzt
 
 type,2
 mat,2
 real,2
 amesh,all
 
 !******************************
 !Simulation definition
 !******************************
 /solu
 /dscale,,1e3
 antype,trans,new
 
 *if,zyx,eq,0,then	!define sinusoidal voltage
  *DIM,spannung,TABLE,1e5,1
 *endif
 *do,i,1,1e5
  spannung(i,0)=tsim/schritte*i
  spannung(i,1)=u0p*SIN(2*pi*f*tsim/schritte*i)
 *enddo
 
 time,tsim
 deltim,tsim/schritte,tsim/schritte,tsim/schritte
 tintp,,0.25,0.5,0.5
 outres,all,all
 
 nsel,s,loc,y,0	!constrain substrate in Y-direction
 nsel,a,loc,y,yhilf+isoueber+laufstr+prdaempf
 d,all,uy,0
 
 lsel,s,,,34,34+(fingernum/2-1)*42,42	!Defining components for plus and minus electrode
 nsll,,1
 cm,plusel,node
 lsel,s,,,55,55+(fingernum/2-1)*42,42
 nsll,,1
 cm,minusel,node
 
 nsel,s,,,plusel	!Applying voltage on electrodes
 d,all,volt,%spannung%
 
 nsel,s,,,minusel
 d,all,volt,0
 
 nsel,s,loc,x,subdick+ged,subdick+ged+id	!Making ground electrode conducting
 cp,1,volt,all
 
 esel,s,type,,1,4,3	!Making steel substrate conducting
 nsle,,1
 cp,2,volt,all
 
 allsel
 parsav,all
 solve	!Solving
 /post1
 !******************************
 !Determining averaged deflection and standard deviation of deflection
 !******************************
 *do,j,1,nint(laufstr/wl-1)	!Determining nodes in undamped traveling region with distance of one wavelength for deflection determination
  nsel,s,loc,y,yhilf+isoueber+j*wl,yhilf+isoueber+(j+1)*wl
  nsel,r,loc,x,subdick
  *get,aktknotpos,node,0,mnloc,y
  nsel,r,loc,y,aktknotpos,aktknotpos+1e-8
  *get,knotnumarr(j),node,0,num,min
 *enddo
 
 set,last
 allsel
 *do,j,1,nint((AnzMW+1)*schritte/tsim/f)	!determining maximum deflection for the <AnzMW> last periods
  *do,k,1,nint(laufstr/wl-1)	!Adressing each node, determined before
   
   *if,j,eq,1,then
    maxux(k)=0
    flagzaehl(k,1)=0	!Don't consider maximum deflection of incomplete first period
    flagzaehl(k,2)=1	!Determining the current period for storing deflection
   *endif
   *if,flagzaehl(k,2),le,AnzMW,then
    *get,aktdis,node,knotnumarr(k),u,x	!Get current deflection
    *if,aktdis,gt,0,then	!Only consider positive deflections
     *if,flagzaehl(k,1),gt,0,then	!Only consider positive deflections of a complete period
      flagzaehl(k,1) = 2
      *if,aktdis,gt,maxux(k),then	!Only consider the maximum
       maxux(k) = aktdis
      *endif
     *endif
    *else
     *if,flagzaehl(k,1),eq,0,then	!Marking the first complete period
      flagzaehl(k,1) = 1
     *elseif,flagzaehl(k,1),eq,2,then	!When it's the first negative deflection, store the maximum deflection and wait for the next positive deflection
      disknot(k,flagzaehl(k,2)) = maxux(k)
      maxux(k)=0
      flagzaehl(k,1) = 1
      flagzaehl(k,2) = flagzaehl(k,2)+1
     *endif
    *endif
   *endif
  *enddo
  set,prev	!select next time step
 *enddo
 
 uxmittel = 0
 uxstdabw = 0
 *do,j,1,AnzMW	!Mittelwert bilden aus allen Verschiebungen	!Determining averaged deflection over all points and considered time steps
  *do,k,1,nint(laufstr/wl-1)
   uxmittel = uxmittel+disknot(k,j)
  *enddo
 *enddo
 uxmittel = uxmittel/AnzMW/nint(laufstr/wl-1)
 
 *do,j,1,AnzMW	!Mittelwert bilden aus allen Verschiebungen	!Determining standard derivation of deflection over all points and considered time steps
  *do,k,1,nint(laufstr/wl-1)
   uxstdabw = uxstdabw+(disknot(k,j)-uxmittel)**2
  *enddo
 *enddo
 uxstdabw = (uxstdabw/AnzMW/nint(laufstr/wl-1))**0.5
 
 !******************************
 !Determining charge and voltage for one period
 !******************************
 fini
 /post1
 *if,zyx,eq,0,then
  *DIM,ladarray,array,AnzMW*schritte/tsim/f,2
 *endif
 *do,g,1,AnzMW*schritte/tsim/f	!consider <AnzMW> periods
  *if,g,eq,1,then
   set,last
  *else
   set,prev
  *endif
  allsel
  nsel,s,,,plusel
  FSUM,0,ALL
  *get,Q_norm,FSUM,0,ITEM,CHRG	!Determining charge on plus electrode
  Ladung = Q_norm*b_IDT
  *get,aktzeit,ACTIVE,0,set,time
  ladarray(g,1) = aktzeit
  ladarray(g,2) = Ladung
 *enddo
 a_l = 0
 *do,g,1,schritte/tsim/f	!Determining first approximation for amplitude of charge
  *if,a_l,lt,ladarray(g,2),then
   a_l=ladarray(g,2)
  *endif
 *enddo
 
 quadfehl = 1e5
 *do,g,1,360
  aktphi_l = g/180*pi
  aktquadfehl=0
  *do,h,1,schritte/tsim/f	!Determining phase shift of charge by mean squared error regarding ideal curve
   aktquadfehl=aktquadfehl+(ladarray(h,2)-a_l*sin(2*pi*f*ladarray(h,1)+aktphi_l))**2
  *enddo
  *if,aktquadfehl,lt,quadfehl,then	!adopt best value
   quadfehl = aktquadfehl
   phi_l = aktphi_l
  *endif
 *enddo
  
 !******************************
 !Fitting amplitude and phase of charge by optimization
 !******************************
 letztquadfehl = 1e5
 letzta_l = a_l
 letztphi_l = phi_l
 schleife = 1
 abbruchzaehler = 0
 variation = 1
 *dowhile,schleife
  *if,abbruchzaehler,eq,20,then	!Reducing variation
   abbruchzaehler = 0
   variation = 0.8*variation
  *endif
  akta_l = (1-RAND(-variation*0.1,variation*0.1))*letzta_l	!random varying of amplitude and phase of charge
  aktphi_l = (1-RAND(-variation*10/180*pi,variation*10/180*pi))*letztphi_l
  
  aktquadfehl=0
  
  *do,h,1,AnzMW*schritte/tsim/f
   aktquadfehl=aktquadfehl+(ladarray(h,2)-akta_l*sin(2*pi*f*ladarray(h,1)+aktphi_l))**2	!Determining mean squared error regarding ideal curve
  *enddo
  
  *if,letztquadfehl,gt,aktquadfehl,then	!adopt best values
   letztquadfehl = aktquadfehl
   letzta_l = akta_l
   letztphi_l = aktphi_l
  *else
   abbruchzaehler = abbruchzaehler+1
  *endif
 
  *if,variation,lt,0.05,then	!Stop, when variation is small enough
   schleife = 0
  *endif
 *enddo
 
 a_l = letzta_l	!finally adopt best values 
 phi_l = letztphi_l
 
 !******************************
 !Calculating impedance
 !******************************
 a_i=a_l*2*pi*f	!Determining current
 phi_i=phi_l+pi/2
 a_z=u0p/a_i	!Determing amplitude and phase of impedance
 phi_z=-phi_i
 R = a_z*cos(phi_z)
 X = a_z*sin(phi_z)
 
 !******************************
 !Storing Results in file
 !******************************
 *cfopen,Deflections_%fingernum%_%wl%_%ctime%,txt,,APPEND
 _f = f
 _fingernum = fingernum
 _mittlux = uxmittel
 _stdabwux = uxstdabw
 _R = R
 _X = X
 *cfwrite,,_f,_fingernum,_mittlux,_stdabwux,_R,_X
 *cfclos
 
 !******************************
 !Store time dependent deflections for each series point
 !******************************
 /POST26 
 FILE,,'rst','.'   
 /UI,COLL,1  
 NUMVAR,200  
 SOLU,191,NCMIT  
 STORE,MERGE 
 FILLDATA,191,,,,1,1 
 REALVAR,191,191
 nsel,s,loc,x,subdick
 nsel,r,loc,y,yhilf+isoueber+wl,yhilf+isoueber+laufstr
 *get,aktknot,NODE,0,num,min
 allsel
 NSOL,2,aktknot,U,X, UX_2   
 STORE,MERGE 
 XVAR,1  
 PLVAR,2,
 *DEL,_P26_EXPORT
 *DIM,_P26_EXPORT,TABLE,schritte,1 
 VGET,_P26_EXPORT(1,0),1 
 VGET,_P26_EXPORT(1,1),2 
 /OUTPUT,'time_dep_displacement_%wl%_%ctime%_%f%','txt','.'   
 *VWRITE,'TIME','UX_2'
 %14C %14C
 *VWRITE,_P26_EXPORT(1,0),_P26_EXPORT(1,1)
 %14.5G %14.5G
 /OUTPUT,TERM
 
 parsav,all
*enddo
save