using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace vcsHDPsurge
{
    public partial class Form1 : Form
    {
        int NN = 30;
        
        public Form1()
        {
            InitializeComponent();
        }
        //---------------------------------------------------------------------------//

        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            main_part();
        }
        //---------------------------------------------------------------------------//
        private void main_part()
        {
            int i,k,ndata,nnn;
            double pp, ph, pw, r, H, z, gamma, kh;
            double lambda, theta, sum, dz;
            System.IO.StreamWriter sw;
            string fnameW="";
            string fname = "";
            string dir="";
            string dat="";
            double[] xx=new double[360];
            double[] yy=new double[360];
            double[] px=new double[360];
            double[] py=new double[360];
            double ratio = 2.0;

            gamma = 1.0;      //水の単位体積重量
            kh = 0.15;        //水平震度
            r = 5.0;          //サージタンク内半径
            H = 10.0;         //水深
            theta = 0.0;      //加速度の方向から半時計回りの角度

            folderBrowserDialog1.Description = "データ出力フォルダを選択してください";
            if (folderBrowserDialog1.ShowDialog() == DialogResult.OK)
            {
                dir = folderBrowserDialog1.SelectedPath;
            }
            
            for (nnn = 0; nnn <= 3; nnn++)
            {
                switch (nnn)
                {
                    case 0: H = 5.0; fname = "hdp_H050.csv"; break;
                    case 1: H = 7.5; fname = "hdp_H075.csv"; break;
                    case 2: H = 10.0; fname = "hdp_H100.csv"; break;
                    case 3: H = 15.0; fname = "hdp_H150.csv"; break;
                }
                fnameW = dir +"\\"+ fname;
                sw=new System.IO.StreamWriter(fnameW,false, System.Text.Encoding.Default);
                dat="#z,pp,ph,pw";
                sw.WriteLine(dat);
                dz = 0.5;
                z=-dz;
                    do{
                        z=z+dz; //サージタンク底面からの距離
                        sum=0.0;
                        for(i=0;i<=NN;i++){
                            lambda=(2.0*(double)i+1.0)/2.0*Math.PI;
                            sum=sum+Math.Pow(-1.0,(double)i)/lambda*FI(r/H,lambda)*Math.Cos(lambda*z/H);
                        }
                        //地震時動水圧強度(速度ポテンシャル)
                        pp=gamma*kh*r*Math.Cos(theta)*sum;
                        //地震時動水圧強度(Housner)
                        ph=Math.Sqrt(3.0)*gamma*kh*H*((H-z)/H-0.5*(H-z)*(H-z)/H/H)*Math.Tanh(Math.Sqrt(3.0)*r/H)*Math.Cos(theta);
                        //地震時動水圧強度(Westergaard)
                        pw=7.0/8.0*gamma*kh*Math.Sqrt(H*(H-z));
                        dat=z.ToString("0.000");
                        dat = dat + "," + pp.ToString("0.000");
                        dat = dat + "," + ph.ToString("0.000");
                        dat = dat + "," + pw.ToString("0.000");
                        sw.WriteLine(dat);
                        if (H - dz < z && z < H + dz) break;
                    }while(z<H);
                sw.Close();
            }

            //動水圧水平分布計算・出力
            H = 15.0; 
            z = 0.0; //サージタンク底面からの距離
            sum=0.0;
            for(i=0;i<=NN;i++){
                lambda=(2.0*(double)i+1.0)/2.0*Math.PI;
                sum=sum+Math.Pow(-1.0,(double)i)/lambda*FI(r/H,lambda)*Math.Cos(lambda*z/H);
            }
            k=0;
            for(i=0;i<=360;i=i+5){
                theta=(double)i/180.0*Math.PI;
                //地震時動水圧強度(速度ポテンシャル)
                pp=gamma*kh*r*Math.Cos(theta)*sum;
                xx[k]=r*Math.Cos(theta);
                yy[k]=r*Math.Sin(theta);
                px[k]=ratio*pp*Math.Cos(theta);
                py[k]=ratio*pp*Math.Sin(theta);
                k=k+1;
            }
            ndata=k-1;

            //動水圧分布出力
            fname = "dtb_H150.csv";
            fnameW = dir + "\\" + fname;
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            dat = "#x,y,dx,dy";
            sw.WriteLine(dat);
            for (k = 0; k <= ndata; k++)
            {
                dat = xx[k].ToString("0.000");
                dat = dat + "," + yy[k].ToString("0.000");
                dat = dat + "," + px[k].ToString("0.000");
                dat = dat + "," + py[k].ToString("0.000");
                sw.WriteLine(dat);
            }
            sw.Close();

            //動水圧成分出力
            fname = "pxy_H150.csv";
            fnameW = dir + "\\" + fname;
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            //x方向（加速度方向）
            dat= "#x-direction";sw.WriteLine(dat);
            dat= "#xx,yy,px,py";sw.WriteLine(dat);
            for (k = 0; k <= ndata; k++)
            {
                if (0.0 <= xx[k])
                {
                    dat = r.ToString("0.000");
                    dat = dat + "," + yy[k].ToString("0.000");
                    dat = dat + "," + px[k].ToString("0.000");
                    dat = dat + "," + (0.0).ToString("0.000");
                    sw.WriteLine(dat);
                }
                if (xx[k] <= 0.0)
                {
                    dat = (-r).ToString("0.000");
                    dat = dat + "," + yy[k].ToString("0.000");
                    dat = dat + "," + px[k].ToString("0.000");
                    dat = dat + "," + (0.0).ToString("0.000");
                    sw.WriteLine(dat);
                }
            }
            //y方向（加速度直交方向）
            dat = "#y-direction"; sw.WriteLine(dat);
            dat = "#xx,yy,px,py"; sw.WriteLine(dat);
            for (k = 0; k <= ndata; k++)
            {
                if (0.0 <= yy[k]) {
                    dat = xx[k].ToString("0.000");
                    dat = dat + "," + r.ToString("0.000");
                    dat = dat + "," + (0.0).ToString("0.000");
                    dat = dat + "," + py[k].ToString("0.000");
                    sw.WriteLine(dat);
                }
                if (yy[k] <= 0.0) { 
                }
                dat = xx[k].ToString("0.000");
                dat = dat + "," + (-r).ToString("0.000");
                dat = dat + "," + (0.0).ToString("0.000");
                dat = dat + "," + py[k].ToString("0.000");
                sw.WriteLine(dat);
            }
            sw.Close();
            
            MessageBox.Show("完了", "完了通知");

        }
        //---------------------------------------------------------------------------//

        double FI(double rbh,double lambda)
        {
            double frbh;
            frbh=2.0/rbh*FMBF(1,lambda*rbh)/(lambda*FMBF(0,lambda*rbh)-FMBF(1,lambda*rbh)/rbh);
            return frbh;
        }
        //---------------------------------------------------------------------------//

        double FMBF(int nu,double z)
        {
            //第1種変形Bessel関数
            double sum;
            int i;
            sum=0.0;
            for(i=0;i<=NN;i++){
                sum=sum+1.0/GAMMAF((double)(i+1))/GAMMAF((double)(nu+i+1))*Math.Pow(0.5*z,(double)(2*i));
            }
            return Math.Pow(0.5*z,(double)nu)*sum;
        }
        //---------------------------------------------------------------------------//

        private double GAMMAF(double z)
        {
            //ガンマ関数
            double a;
            a=z/Math.E*Math.Sqrt(z*Math.Sinh(1.0/z)+1.0/810.0/Math.Pow(z,6.0));
            return Math.Sqrt(2.0*Math.PI/z)*Math.Pow(a,z);
        }
        //---------------------------------------------------------------------------//
    
    }
}