/*******************************************/
/* Vibmdof */
/* All Rights Reserved by s.fukuda */
/* From the equation */
/* [ A + jB ]{ x }={ P } */
/* */
/* the solution { x } may be written */
/* { x }=[ C + jD ]{ P } */
/* */
/* =[ E + jF ] */
/*******************************************/
import java.applet.Applet;
import java.awt.event.*;
import java.awt.*;
import java.lang.*;
import java.util.*;

public class vibmdof1 extends Applet {
String msg;
Graphics g;

/* Edit following nm,ns,nd to change numbers of input panel */
int nm=10;
int ns=10;
int nd=10;
int idispm;

TextField nnm,nns,nnd,dispm;
String xnnm,xnns,xnnd,xdispm;

TextField mas[]=new TextField[nm];
TextField foc[]=new TextField[nm];
TextField stf[]=new TextField[ns];
TextField dmp[]=new TextField[nd];
TextField stf_e1[]=new TextField[ns];
TextField stf_e2[]=new TextField[ns];

int xxm[]=new int[nm];
int xfc[]=new int[nm];
int xxs[]=new int[ns];
int xdm[]=new int[nd];
int xxs_e1[]=new int[ns];
int xxs_e2[]=new int[ns];

Thread th=null;

public void init(){
msg=getParameter("message");
g=getGraphics();

Panel p=new Panel();
Panel p0=new Panel();
Panel p1=new Panel();
Panel p2=new Panel();
Panel p3=new Panel();
Panel p4=new Panel();
Panel p5=new Panel();
Panel p6=new Panel();
Panel p10=new Panel();
Panel p20=new Panel();
Panel p50=new Panel();
Panel p60=new Panel();
p.setLayout(new GridLayout(5,1));
p0.setLayout(new GridLayout(5,1));
p1.setLayout(new GridLayout(nm,1));
p2.setLayout(new GridLayout(ns,1));
p3.setLayout(new GridLayout(ns,1));
p4.setLayout(new GridLayout(ns,1));
p5.setLayout(new GridLayout(nd,1));
p6.setLayout(new GridLayout(nm,1));
p10.setLayout(new GridLayout(2,1));
p20.setLayout(new GridLayout(2,1));
p50.setLayout(new GridLayout(2,1));
p60.setLayout(new GridLayout(nm,1));


p0.add(new Button("Calc"));
p0.add("",nnm=new TextField("1",2));
p0.add("",nns=new TextField("1",2));
p0.add("",nnd=new TextField("1",2));
p0.add("",dispm=new TextField("0",2));

p1.add("",mas[0]=new TextField("10",3));
for (i=1; i<nm; i++){
p1.add("",mas[i]=new TextField("0",3));
}
p2.add("",stf[0]=new TextField("10000000",6));
p3.add("",stf_e1[0]=new TextField(Integer.toString(0),2));
p4.add("",stf_e2[0]=new TextField(Integer.toString(-1),2));

for (i=1; i<ns; i++){
p2.add("",stf[i]=new TextField("0",6));
p3.add("",stf_e1[i]=new TextField(Integer.toString(0),2));
p4.add("",stf_e2[i]=new TextField(Integer.toString(0),2));
}

p5.add("",dmp[0]=new TextField("100",3));
for (i=1; i<nd; i++){
p5.add("",dmp[i]=new TextField("0",3));
}
p6.add("",foc[0]=new TextField("1000",3));
for (i=1; i<nm; i++){
p6.add("",foc[i]=new TextField("0",3));
}

p.add("Center",new Label("Click >>"));
p.add("Center",new Label("質量の数"));
p.add("Center",new Label("ばねの数"));
p.add("Center",new Label("減衰の数"));
p.add("Center",new Label("表示する質量"));

p10.add("Center",new Label(" Mass"));
p10.add("Center",new Label(" (kg)"));
p20.add("Center",new Label(" Stiff."));
p20.add("Center",new Label(" (N/m)"));
p50.add("Center",new Label(" Damp."));
p50.add("Center",new Label("(N*S/m)"));
p60.add("Center",new Label(" Force"));
p60.add("Center",new Label(" ( N )"));

add("Center",p);
add("Center",p0);

add("Center",p10);
add("Center",p1);
add("Center",p60);
add("Center",p6);
add("Center",p20);
add("Center",p2);
add("Center",p50);
add("Center",p5);

add("Center",new Label(" Conect to"));
add("Center",p3);
add("Center",p4);
}
/**
public void actionPerformed(ActionEvent e){
if ("Calc".equals(e.getActionCommand())){
**/
public boolean action(Event evt, Object arg) {
if (evt.target instanceof Button) {
repaint();
}
return true;
}

/* dimension */
double mass[]=new double[20];
double stiff[][]=new double[20][3];
double force[]=new double[20];
double damp[][]=new double[20][3];
double k0[][]=new double[20][20];
double c0[][]=new double[20][20];
double m1[][]=new double[40][80];
double m2[][]=new double[40][40];
double e[]=new double[20];
double f[]=new double[20];
double m3[][]=new double[40][1024];
double m4[][]=new double[20][1024];
double m5[][]=new double[20][1024];
double i0[][]=new double[40][20];

double pi=3.14159;
double d1,s1,w,pivi,w1,w3,t;
double dmax,dmin,f1,f2;
double xscale,yscale,ds,ss;
double ay1,ay2,aywind,pywind,py1,py2,x1,x2;
int i,j,k,k1,nj,e1,e2,freqpnt,fp;
int xwind,ywind;
int xmin,xmax,aymax,aymin,pymax,py0,pymin;
int ii,jj;

public void paint(Graphics g){

g.setColor(Color.gray);
g.fillRect(0,250,760,600);
g.setColor(Color.black);
g.fillRect(440,300,310,200);

String xm[]=new String[20];
String xs[][]=new String[20][3];
String xd[]=new String[20];
String xf[]=new String[20];

xnnm=nnm.getText();
xnns=nns.getText();
xnnd=nnd.getText();
xdispm=dispm.getText();

nm=Integer.parseInt((String)xnnm);
ns=Integer.parseInt((String)xnns);
nd=Integer.parseInt((String)xnnd);
idispm=Integer.parseInt((String)xdispm);

for (i=0; i<nm; i++){
xm[i]=mas[i].getText();
}

for (i=0; i<ns; i++){
xs[i][0]=stf[i].getText();
xs[i][1]=stf_e1[i].getText();
xs[i][2]=stf_e2[i].getText();
}

for (i=0; i<nd; i++){
xd[i]=dmp[i].getText();
}

for (i=0; i<nm; i++){
xf[i]=foc[i].getText();
}

for (i=0; i<nm; i++){
xxm[i]=Integer.parseInt((String)xm[i]);
}

for (i=0; i<ns; i++){
xxs[i]=Integer.parseInt((String)xs[i][0]);
xxs_e1[i]=Integer.parseInt((String)xs[i][1]);
xxs_e2[i]=Integer.parseInt((String)xs[i][2]);
}

for (i=0; i<nd; i++){
xdm[i]=Integer.parseInt((String)xd[i]);
}

for (i=0; i<nm; i++){
xfc[i]=Integer.parseInt((String)xf[i]);
}

Color col1=Color.white;
Color col2=Color.blue;
Color col3=Color.yellow;
Color col4=Color.red;
g.setColor(col1);
/*
ywind=(int)size().height;
xwind=(int)size().width;
*/
xwind=640;
ywind=800;

fp=256;
f1=1;
f2=1000;

for (i=0; i<nm; i++){
mass[i]=(double)xxm[i];
}

for (i=0; i<ns; i++){
stiff[i][0]=(double)xxs[i];
stiff[i][1]=(double)xxs_e1[i];
stiff[i][2]=(double)xxs_e2[i];
}

for (i=0; i<ns; i++){
damp[i][0]=(double)xdm[i];
damp[i][1]=(double)xxs_e1[i];
damp[i][2]=(double)xxs_e2[i];
}

for (i=0; i<nm; i++){
force[i]=(double)xfc[i];
}

g.setColor(Color.red);
g.drawString("VibMdof (v1.3.0)",10,15);
g.drawString(" All rights reserved.",20,30);
g.setColor(Color.blue);
g.drawString("パラメータを変更し",10,45);
g.drawString("Calc をクリックして再計算",20,60);

g.setColor(Color.blue);
g.drawString("表示する質量番号を選択し",20,200);
g.drawString("その質量の周波数応答の",20,215);
g.drawString("グラフを表示します。",20,230);
g.drawString("( 0 =< M# < 質量の数 )",20,245);
/*
g.drawString("Mass",550,270);
for (i=0; i<nm; i++){
g.drawString("m"+i+" = "+Double.toString(mass[i]),550,15*(i+1)+270);
}
g.drawString("Stiffness ( m-1 means fixed end ) ",550,270+15*(nm+2));
for (i=0; i<ns; i++){
g.drawString(Double.toString(stiff[i][0])+"(conected to m"+Double.toString(stiff[i][1])+" - m"+Double.toString(stiff[i][2])+")",550,15*(i+3+nm)+270);
}
g.drawString("Damping",550,270+15*(nm+ns+4));
for (i=0; i<nd; i++){
g.drawString(Double.toString(damp[i][0])+"(conected to m"+Double.toString(damp[i][1])+" - m"+Double.toString(damp[i][2])+")",550,15*(i+5+nm+ns)+270);
}
g.drawString("Force",550,270+15*(nm+ns+nd+6));
for (i=0; i<nm; i++){
g.drawString("force(m"+i+")="+Double.toString(force[i]),550,15*(i+9+2*nm+ns)+270);
}
*/


/* Create system matrix [ K ] */
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
k0[i][j]=0.0;
}
}
for (i=0; i<ns; i++){
e1=(int)(stiff[i][1]);
e2=(int)(stiff[i][2]);
if (e1==-1)
k0[e2][e2]=k0[e2][e2]+stiff[i][0];
else if (e2==-1)
k0[e1][e1]=k0[e1][e1]+stiff[i][0];
else{
k0[e1][e2]=k0[e1][e2]-stiff[i][0];
k0[e2][e1]=k0[e2][e1]-stiff[i][0];
k0[e1][e1]=k0[e1][e1]+stiff[i][0];
k0[e2][e2]=k0[e2][e2]+stiff[i][0];
}
}
/* printing */
/*
g.setColor(col2);
g.drawString("Stiffness Matrix",270,50);
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
g.drawString(Double.toString(k0[i][j]),100*j+270,15*(i+1)+50);
}
}
*/

/* Create system matrix [ C ] */
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
c0[i][j]=0.0;
}
}
for (i=0; i<ns; i++){
e1=(int)(damp[i][1]);
e2=(int)(damp[i][2]);
if (e1==-1)
c0[e2][e2]=c0[e2][e2]+damp[i][0];
else if (e2==-1)
c0[e1][e1]=c0[e1][e1]+damp[i][0];
else{
c0[e1][e2]=c0[e1][e2]-damp[i][0];
c0[e2][e1]=c0[e2][e1]-damp[i][0];
c0[e1][e1]=c0[e1][e1]+damp[i][0];
c0[e2][e2]=c0[e2][e2]+damp[i][0];
}
}

/* printing */
/*
g.setColor(col2);
g.drawString("Damping Matrix",100,15*(nm+1)+50);
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
g.drawString(Double.toString(c0[i][j]),100*j+100,15*(i+nm+2)+50);
}
}
*/

/****************/
/* main routine */
/****************/
w1=2*pi*f1;
s1=w1;
d1=(f2-f1)/(fp-1);
w3=2*pi*d1;
for (j=0; j<nm; j++){
for (i=0; i<2*nm; i++){
i0[ i ][ j ]=0;

i0[ j ][ j ]=1;
}
}


/* frequency response function */
for (freqpnt=0; freqpnt<fp; freqpnt++){
w1=s1+(freqpnt)*w3;

/* Create system matrix [ M1 ]=[ A ] */
/* [ B ] */
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
if (i==j)
m1[i][i]=k0[i][i]-w1*w1*mass[i]; /* system matrix [ A ] */
else{
m1[i][j]=k0[i][j];
m1[i][i]=k0[i][i]-w1*w1*mass[i]; /* system matrix [ A ] */
}
}
for (k=0; k<nm; k++){
m1[nm+i][k]=w1*c0[i][k]; /* system matrix [ B ] */
}
}

/* Inversion of Complex matrix */
/* sweep out method [ A -B I 0 ] */
/* [ B A 0 I ] */

for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
m1[i][nm+j]= -m1[nm+i][j];
m1[nm+i][nm+j]= m1[i][j];
}
}
for (i=0; i<2*nm; i++){
nj= 2*nm;
for (j=nj; j<4*nm; j++){
m1[i][j]=0;
}
m1[i][2*nm+i]=1;
}
for (k=0; k<2*nm; k++){
pivi=m1[k][k];
k1=k+1;
for (j=k1; j<4*nm; j++){
m1[k][j]=(m1[k][j])/pivi;
}
for (i=0; i<2*nm; i++){
if(i !=k){
w=m1[i][k];
for (j=k1; j<4*nm; j++){
m1[i][j]=m1[i][j]-w*m1[k][j];
}
}
}
}
/* end of inversion */


/* calculate M2 */
for (k=0; k<nm; k++){
for (i=0; i<2*nm; i++){
t=0;
for (j=0; j<2*nm; j++){
t=t+m1[i][2*nm+j]*i0[j][k];
}
m2[i][k]=t;
}
}
for (i=0; i<nm; i++){
t=0;
for (j=0; j<nm; j++){
t=t+m2[i][j]*force[j];
}
e[i]=t;
}
for (i=0; i<nm; i++){
t=0;
for (j=0; j<nm; j++){

t=t+m2[nm+i][j]*force[j];
}
f[i]=t;
}
for (i=0; i<nm; i++){
m3[i][freqpnt]=e[i]*w1*w1;
m3[nm+i][freqpnt]=f[i]*w1*w1;
}
} /* end of frequency response function */


/* bode plot */
/* amplitude */
for (i=0; i<nm; i++){
for (j=0; j<fp; j++){
m4[i][j]=Math.sqrt((m3[i][j])*(m3[i][j])+(m3[nm+i][j])*(m3[nm+i][j]));
}
}
/*
for (i=0; i<nm; i++){
g.drawString("m"+i ,470+100*i,100);
for (j=0; j<fp; j++){
g.drawString(Double.toString(m4[i][j]),500+100*i,j*15+100);
}
}
*/
/* phase */
for (i=0; i<nm; i++){
for (j=0; j<fp; j++){
m5[i][j]=180/pi*Math.atan2( m3[nm+i][j],m3[i][j] );
}
}

/* drawing routine */

i=idispm;

dmax=1000;
dmin=1;

d1=(f2-f1)/(fp-1);
xmin=60;
xmax=400;
/**
aymax=390+300*i;
aymin=510+300*i;
pymax=270+300*i;
py0= 300+300*i;
pymin=330+300*i;
**/
aymax=360;
aymin=480;
pymax=270;
py0= 300;
pymin=330;

g.setColor(Color.white);
g.drawString("Phase",xmin+10,pymax-5);
g.setColor(Color.yellow);
g.drawString("M"+i+" is selected",xmin+80,pymax-5);

g.setColor(col1);
g.drawLine(xmin,aymin, xmax,aymin );
g.drawLine(xmin,aymin, xmin,aymax );
g.drawLine(xmin,py0, xmax,py0 );
g.drawLine(xmin,pymin, xmin,pymax );
g.drawString("Amplitude",xmin+10,aymax-10);
g.drawString("Impulse response",xmax+50,pymax+10);
g.drawString("Frequency (Hz)",xmin+50,aymin+20);
g.drawString(Double.toString(f1),xmin-10,aymin+20);
g.drawString(Double.toString(f2),xmax-10,aymin+20);

g.drawString(Double.toString(dmin),xmin-40,aymin+10);
g.drawString(Double.toString(dmax),xmin-40,aymax+10);
g.drawString("Acc.",xmin-40,(aymin+aymax)/2);
g.drawString("(m/s^2)",xmin-50,(aymin+aymax)/2+20);
// g.drawString("Phase",xmax+10,py0+10);
g.drawString(" 0",xmin-30,py0+10);
g.drawString(" 180",xmin-30,pymax+10);
g.drawString("-180",xmin-30,pymin+10);

/* scale of amplitude */
xwind=xmax-xmin;
ywind=aymin-aymax;
for (ds=1; ds<=8; ds++){
for (ss=1; ss<=10; ss++){
yscale=aymin-ywind*Math.log(Math.pow(10,(ds-1))*ss/dmin)/Math.log(dmax/dmin);

if(yscale <= aymin && yscale >=aymax ){
g.drawLine((int)xmin,(int)yscale,(int)xmax,(int)yscale );
}
}
}

for (ds=1; ds<=8; ds++){
for (ss=1; ss<=10; ss++){
xscale=xmin+xwind*Math.log(Math.pow(10,(ds-1))*ss/f1)/Math.log(f2/f1);
if(xscale >= xmin && xscale <=xmax ){
g.drawLine((int)xscale,(int)aymin, (int)xscale,(int)aymax );
g.drawLine((int)xscale,(int)py0, (int)xscale,(int)(py0-5) );
}
}
}


/* amplitude */
g.setColor(Color.blue);
x1 =xmin;
if( m4[i][0]==0){
ay1=aymin;
}
else {
ay1=aymin-ywind*Math.log (( m4[i][0]/dmin))/Math.log(dmax/dmin);
if (ay1 > aymin){
ay1=aymin;
g.setColor(Color.white);
}
}

for (j=1; j<fp; j++){
ay2=aymin-ywind*Math.log(( m4[i][j]/dmin))/Math.log(dmax/dmin);
x2 =xmin+xwind*Math.log((f1+(j)*d1)/f1)/Math.log(f2/f1);
if (ay2 > aymin){
ay2=aymin;
g.setColor(Color.white);
}
else{
g.setColor(Color.blue);
}
g.drawLine((int)x1,(int)ay1,(int)x2,(int)ay2 );
x1=x2;
ay1=ay2;
}

/* phase */
g.setColor(col3);
pywind=py0-pymax;
py1=py0-pywind*m5[i][0]/180;
x1=xmin;

for (j=1; j<fp; j++){
py2=py0-pywind*m5[i][j]/180;
x2=xmin+xwind*Math.log((f1+(j)*d1)/f1)/Math.log(f2/f1);

g.drawLine((int)x1,(int)py1,(int)x2,(int)py2 );
x1=x2;
py1=py2;
}


/************************************************/
/* Calculate Impulse response using inverse FFT.*/
/* if f=-1.0 then fft if f= 1.0 then ifft */
/************************************************/
double f;
int i,i0,i1,j,l1,n,ns,n1,arg,arg2;
double x[]=new double[1024];
double y[]=new double[1024];
int m[]=new int[1024];
int gx0,gx1,gx2,gy0,gy1,gy2,dyw,xwind,ywind,xposi,yposi;
double dy1,dy2,ymax,ydam,amp,phase;
double s,c,sc,x1,y1,t;
g.setColor(col1);
xwind=290;
ywind=600;
xposi=450;
yposi=250;

/* input data */
/* note: change dimension of arrays */
n=512;
i=idispm;
for (j=0; j<n/2; j++){
x[n/2+j]=m3[i][j];
y[n/2+j]=m3[nm+i][j];
}

/* Hermitian's symmetricity */
/*
for (j=1; j<n/2; j++){
x[n/2 -j]= x[n/2 +j];
y[n/2 -j]=-y[n/2 +j];
}
x[0]=0.0;
y[0]=0.0;
x[n/2]=0.0;
y[n/2]=0.0;
for (j=1; j<n; j++){
y[j]=- y[j];
}
*/
/*******/
/* fft */
/*******/
/* if f=-1.0 then fft if f= 1.0 then rfft */
f =1.0;
n1=n/2;
ns=n/2;
sc=2*3.14159/(double)n;

while (ns>=1) /* main loop */
{
for (l1=0; l1<n; l1+=(2*ns))
{
arg=m[l1]/2;
arg2=arg+n1;
c=Math.cos(sc*(double)arg);
s=Math.sin(f*sc*(double)arg);
for (i0=l1; i0<l1+ns; i0++)
{
i1 =i0+ns;
x1 =x[i1]*c-y[i1]*s; y1 =y[i1]*c+x[i1]*s;
x[i1]=x[i0]-x1; y[i1]=y[i0]-y1;
x[i0]=x[i0]+x1; y[i0]=y[i0]+y1;
m[i0]=arg; m[i1]=arg2;
}
}
ns=ns/2;
}

if (f<0.0){
for (i=0; i<n; i++){
x[i]/=(double)n;
y[i]/=(double)n;
}
}

for (i=0; i<n; i++) { /* bit reverse operation */
if ( (j=m[i]) > i ){
t=x[i]; x[i]=x[j]; x[j]=t;
t=y[i]; y[i]=y[j]; y[j]=t;
}
}

/* draw output data */
g.setColor(Color.blue);
g.drawLine(xposi,ywind/4 + yposi,xwind+xposi,ywind/4 + yposi);

g.setColor(Color.green);
ymax=x[0];
for (i=0; i<n; i++){
if( x[i]<0.0 ){
ydam=-x[i];
}
else {
ydam= x[i];
}
if( ymax < ydam ){
ymax=ydam;
}
}
gy0=ywind/4 + yposi;
gx1=xposi;
dy1=x[0]/ymax;
dyw=ywind/8;
gy1=-(int)(dy1 * dyw) +gy0;
for (i=1; i<n; i++){
gx2= i * (xwind)/(n-1)+xposi;
dy2=x[i]/ymax;
gy2=-(int)(dy2 * dyw) +gy0; /* real only */
g.drawLine( gx1,gy1,gx2,gy2);
gx1=gx2;
gy1=gy2;
}

g.setColor(Color.white);
showStatus(msg);
} /* end of paint */
} /* end of class vibmdof */