------------------------data1----------------------------
ゲイン gain 位相 phase
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 0 次数 0 係数（次数が低い順） 1.0
分母の表現方法 0 次数 3 係数（次数が低い順） 1.0 11.1 11.1 1.0
------------------------data2----------------------------
ゲイン gain 位相 phase
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 1 次数 0 係数（次数が低い順） 1.0
分母の表現方法 1 次数 3 係数（次数が低い順） 1.0 0.1 1.0 1.0 1.0 10.0
------------------------プログラム-----------------------
/********************************/
/* 伝達関数のゲインと位相の計算 */
/* coded by Y.Suganuma */
/********************************/
import java.io.*;
import java.util.*;
public class Test {
/****************/
/* main program */
/****************/
public static void main(String args[]) throws IOException
{
double f, ff, fl, fu, h, gain, phase = 0.0, x, uc;
double bo[] = new double [1];
double si[] = new double [1];
int i1, k, k1, ms, mb, m, n;
Complex g;
String str;
StringTokenizer token;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
uc = 90.0 / Math.asin(1.0);
/*
データの入力
*/
// 出力ファイル名とファイルのオープン
str = in.readLine();
token = new StringTokenizer(str, " ");
token.nextToken();
PrintStream g_o = new PrintStream(new FileOutputStream(token.nextToken()));
token.nextToken();
PrintStream p_o = new PrintStream(new FileOutputStream(token.nextToken()));
// 伝達関数データ
str = in.readLine();
token = new StringTokenizer(str, " ");
token.nextToken();
fl = Double.parseDouble(token.nextToken());
fu = Double.parseDouble(token.nextToken());
token.nextToken();
k = Integer.parseInt(token.nextToken());
// 分子
str = in.readLine();
token = new StringTokenizer(str, " ");
token.nextToken();
ms = Integer.parseInt(token.nextToken());
token.nextToken();
m = Integer.parseInt(token.nextToken());
token.nextToken();
// 因数分解されていない
if (ms == 0) {
k1 = m + 1;
si = new double [k1];
for (i1 = 0; i1 < k1; i1++)
si[i1] = Double.parseDouble(token.nextToken());
}
// 因数分解されている
else {
k1 = (m == 0) ? 1 : 2 * m;
si = new double [k1];
for (i1 = 0; i1 < k1; i1++)
si[i1] = Double.parseDouble(token.nextToken());
}
// 分母
str = in.readLine();
token = new StringTokenizer(str, " ");
token.nextToken();
mb = Integer.parseInt(token.nextToken());
token.nextToken();
n = Integer.parseInt(token.nextToken());
token.nextToken();
// 因数分解されていない
if (mb == 0) {
k1 = n + 1;
bo = new double [k1];
for (i1 = 0; i1 < k1; i1++)
bo[i1] = Double.parseDouble(token.nextToken());
}
// 因数分解されている
else {
k1 = (n == 0) ? 1 : 2 * n;
bo = new double [k1];
for (i1 = 0; i1 < k1; i1++)
bo[i1] = Double.parseDouble(token.nextToken());
}
/*
ゲインと位相の計算
*/
h = (log10(fu) - log10(fl)) / k;
ff = log10(fl);
for (i1 = 0; i1 <= k; i1++) {
// 周波数の対数を元に戻す
f = Math.pow(10.0, ff);
// 値の計算
g = G_s(f, ms, m, si, mb, n, bo);
// ゲインと位相の計算
gain = 20.0 * log10(Complex.abs(g));
x = Complex.angle(g) * uc;
while (Math.abs(phase-x) > 180.0) {
if (x-phase > 180.0)
x -= 360.0;
else {
if (x-phase < -180.0)
x += 360.0;
}
}
phase = x;
// 出力
g_o.println(f + " " + gain);
p_o.println(f + " " + phase);
// 次の周波数
ff += h;
}
g_o.close();
p_o.close();
}
/************/
/* log10(x) */
/************/
static double log10(double x)
{
return Math.log(x) / Math.log(10.0);
}
/********************************************************************/
/* 伝達関数のｓにｊωを代入した値の計算 */
/* ff : ω（周波数） */
/* ms : 分子の表現方法 */
/* =0 : 多項式（si[0]+si[1]*s+si[1]*s^2+･･･） */
/* =1 : 因数分解（(si[0]+si[1]*s)*(si[2]+si[3]*s)*･･･） */
/* m : 分子の次数 */
/* si : 分子多項式の係数 */
/* mb : 分母の表現方法 */
/* =0 : 多項式（bo[0]+bo[1]*s+bo[1]*s^2+･･･） */
/* =1 : 因数分解（(bo[0]+bo[1]*s)*(bo[2]+bo[3]*s)*･･･） */
/* n : 分母の次数 */
/* bo : 分母多項式の係数 */
/* return : 結果 */
/********************************************************************/
static Complex G_s(double ff, int ms, int m, double si[], int mb, int n, double bo[])
{
Complex f, x, y;
// 周波数を複素数に変換
f = new Complex (0.0, ff);
// 分子
x = value(f, ms, m, si);
// 分母
y = value(f, mb, n, bo);
return Complex.dev(x, y);
}
/****************************************************************/
/* 多項式のｓにｊωを代入した値の計算 */
/* f : ｊω（周波数，複素数） */
/* ms : 多項式の表現方法 */
/* =0 : 多項式（a[0]+a[1]*s+a[1]*s^2+･･･） */
/* =1 : 因数分解（(a[0]+a[1]*s)*(a[2]+a[3]*s)*･･･） */
/* n : 多項式の次数 */
/* a : 多項式の係数 */
/* return : 結果 */
/****************************************************************/
static Complex value(Complex f, int ms, int n, double a[])
{
int i1, k1;
Complex x;
// 因数分解されていない
if (ms == 0) {
x = new Complex (0.0, 0.0);
for (i1 = 0; i1 <= n; i1++)
x = Complex.add(x, Complex.mul(new Complex(a[i1]), Complex.pow(f, i1)));
}
// 因数分解されている
else {
if (n == 0)
x = new Complex (a[0]);
else {
x = new Complex (1.0, 0.0);
k1 = 2 * n;
for (i1 = 0; i1 < k1; i1 += 2)
x = Complex.mul(x, Complex.add(new Complex(a[i1]), Complex.mul(new Complex(a[i1+1]), f)));
}
}
return x;
}
}
/****************************/
/* クラスComplex */
/* coded by Y.Suganuma */
/****************************/
class Complex {
private double r; // 実部
private double i; // 虚部
/******************/
/* コンストラクタ */
/* a : 実部 */
/* b : 虚部 */
/******************/
Complex(double a, double b)
{
r = a;
i = b;
}
/******************/
/* コンストラクタ */
/* a : 実部 */
/******************/
Complex(double a)
{
r = a;
i = 0.0;
}
/******************/
/* コンストラクタ */
/******************/
Complex()
{
r = 0.0;
i = 0.0;
}
/********/
/* 出力 */
/********/
void output(PrintStream out, Complex a)
{
out.print("(" + a.r + ", " + a.i + ")");
}
/******************/
/* 複素数の絶対値 */
/******************/
static double abs(Complex a)
{
double x;
x = Math.sqrt(a.r * a.r + a.i * a.i);
return x;
}
/****************/
/* 複素数の角度 */
/****************/
static double angle(Complex a)
{
double x;
if (a.r == 0.0 && a.i == 0.0)
x = 0.0;
else
x = Math.atan2(a.i, a.r);
return x;
}
/****************/
/* 複素数のｎ乗 */
/****************/
static Complex pow(Complex a, int n)
{
int i1;
Complex c = new Complex (1);
if (n >= 0) {
for (i1 = 0; i1 < n; i1++)
c = Complex.mul(c, a);
}
else {
for (i1 = 0; i1 < -n; i1++)
c = Complex.mul(c, a);
c = Complex.dev(new Complex(1.0), c);
}
return c;
}
/****************/
/* 複素数の加算 */
/****************/
static Complex add(Complex a, Complex b)
{
Complex c = new Complex();
c.r = a.r + b.r;
c.i = a.i + b.i;
return c;
}
/****************/
/* 複素数の減算 */
/****************/
static Complex sub(Complex a, Complex b)
{
Complex c = new Complex();
c.r = a.r - b.r;
c.i = a.i - b.i;
return c;
}
/****************/
/* 複素数の乗算 */
/****************/
static Complex mul(Complex a, Complex b)
{
Complex c = new Complex();
c.r = a.r * b.r - a.i * b.i;
c.i = a.i * b.r + a.r * b.i;
return c;
}
/****************/
/* 複素数の除算 */
/****************/
static Complex dev(Complex a, Complex b)
{
double x;
Complex c = new Complex();
x = 1.0 / (b.r * b.r + b.i * b.i);
c.r = (a.r * b.r + a.i * b.i) * x;
c.i = (a.i * b.r - a.r * b.i) * x;
return c;
}
}