Option Explicit On
Option Strict On
Public Class Form1

    Private NDmax As Integer = 100
    Private strcom As String       ' Comment 
    Private nd As Integer              ' Numbers of calculation 
    Private IE_inp(NDmax) As Integer ' 0:Inner pressure, 1:External pressure 
    Private PP_inp(NDmax) As Double ' Water pressure 
    Private TT_inp(NDmax) As Double ' Temperature change (+:increse of temperature) 
    Private aa_inp(NDmax) As Double ' Internal radius of pressure tunnel 
    Private bb_inp(NDmax) As Double ' Excavation radius 
    Private rr_inp(NDmax) As Double ' Model edge radius 
    Private cc_inp(NDmax) As Double ' Cover for reinforcement 
    Private ta_inp(NDmax) As Double ' Equivalent thickness of inner reinforcement 
    Private tb_inp(NDmax) As Double ' Equivalent thickness of outer reinforcement 
    Private Ec_inp(NDmax) As Double ' Elastic modulus of concrete 
    Private nc_inp(NDmax) As Double ' Poisson ratio of concrete 
    Private ac_inp(NDmax) As Double ' Thermal expansion coefficient of concrete 
    Private Es_inp(NDmax) As Double ' Elastic modulus of reinforcement 
    Private ns_inp(NDmax) As Double ' Poisson ratio of reinforcement 
    Private as_inp(NDmax) As Double ' Thermal expansion coefficient of reinforcement 
    Private Eg_inp(NDmax) As Double ' Elastic modulus of bed rock 
    Private ng_inp(NDmax) As Double ' Poisson ratio of bed rock 

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim fnameR As String = ""
        Dim fnameW As String = ""
        Dim dat As String
        Dim sbuf() As String
        Dim delim As Char = ","c
        Dim k As Integer
        Dim nd As Integer

        OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory()
        If OpenFileDialog1.ShowDialog() = DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        If SaveFileDialog1.ShowDialog() = DialogResult.OK Then fnameW = SaveFileDialog1.FileName

        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        strcom = sr.ReadLine()
        k = 0
        Do Until sr.EndOfStream
            dat = sr.ReadLine()
            If dat.Substring(0, 1) <> "#" Then
                sbuf = dat.Split(delim)
                IE_inp(k) = Integer.Parse(sbuf(0))
                PP_inp(k) = Double.Parse(sbuf(1))
                TT_inp(k) = Double.Parse(sbuf(2))
                aa_inp(k) = Double.Parse(sbuf(3))
                bb_inp(k) = Double.Parse(sbuf(4))
                rr_inp(k) = Double.Parse(sbuf(5))
                cc_inp(k) = Double.Parse(sbuf(6))
                ta_inp(k) = Double.Parse(sbuf(7))
                tb_inp(k) = Double.Parse(sbuf(8))
                Ec_inp(k) = Double.Parse(sbuf(9))
                nc_inp(k) = Double.Parse(sbuf(10))
                ac_inp(k) = Double.Parse(sbuf(11))
                Es_inp(k) = Double.Parse(sbuf(12))
                ns_inp(k) = Double.Parse(sbuf(13))
                as_inp(k) = Double.Parse(sbuf(14))
                Eg_inp(k) = Double.Parse(sbuf(15))
                ng_inp(k) = Double.Parse(sbuf(16))
                k = k + 1
            End If
        Loop
        sr.Close()
        nd = k - 1

        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        sw.WriteLine(strcom)
        sw.WriteLine("*Input data")
        sw.WriteLine("k,IE,PP,TT,aa,bb,rr,cc,ta,tb,Ec,nc,ac,Es,ns,as,Eg,ng")
        For k = 0 To nd
            dat = IE_inp(k).ToString()
            dat = dat & "," & PP_inp(k).ToString("#.###")
            dat = dat & "," & TT_inp(k).ToString("#.###")
            dat = dat & "," & aa_inp(k).ToString("#.###")
            dat = dat & "," & bb_inp(k).ToString("#.###")
            dat = dat & "," & rr_inp(k).ToString("#.###")
            dat = dat & "," & cc_inp(k).ToString("#.###")
            dat = dat & "," & ta_inp(k).ToString("#.###")
            dat = dat & "," & tb_inp(k).ToString("#.###")
            dat = dat & "," & Ec_inp(k).ToString("#.###")
            dat = dat & "," & nc_inp(k).ToString("#.###")
            dat = dat & "," & ac_inp(k).ToString("#.###")
            dat = dat & "," & Es_inp(k).ToString("#.###")
            dat = dat & "," & ns_inp(k).ToString("#.###")
            dat = dat & "," & as_inp(k).ToString("#.###")
            dat = dat & "," & Eg_inp(k).ToString("#.###")
            dat = dat & "," & ng_inp(k).ToString("#.###")
            sw.WriteLine(dat)
        Next k

        sw.WriteLine("*Output data\n")
        ' k     : Sequence                                                     
        ' IE    : Inner pressure (=0) or Outer pressure (=1)                   
        ' Eg    : Water pressure                                               
        ' dT    : Temperature change                                           
        ' sr_c  : Stress in the r-direction at inner surface concrete             
        ' st_c  : Stress in the t-direction at inner surface concrete             
        ' sr_si : Stress in the r-direction of inner reinforcement                
        ' st_si : Stress in the t-direction of inner reinforcement                
        ' sr_so : Stress in the r-direction of outer reinforcement                
        ' st_so : Stress in the t-direction of outer reinforcement                
        ' sr_g  : Stress in the r-direction at outer surface concrete or bed rock 
        ' st_g  : Stress in the t-direction at outer surface concrete or bed rock 
        ' ua    : Radial displacement of inner surface concrete                
        ' ub    : Radial displacement of outer surface concrete                
        sw.WriteLine("k,IE,Eg,dT,sr_c,st_c,sr_si,st_si,sr_so,st_so,sr_g,st_g,ua,ub")

        ' Execution of calculation 
        For k = 0 To nd
            Select Case IE_inp(k)
                Case 0
                    If tb_inp(k) <= 0 Then Call main_sng_pin(sw, k)
                    If 0 < tb_inp(k) Then Call main_dbl_pin(sw, k)
                Case 1
                    If tb_inp(k) <= 0 Then Call main_sng_pex(sw, k)
                    If 0 < tb_inp(k) Then Call main_dbl_pex(sw, k)
            End Select
        Next k

        sw.Close()

        MessageBox.Show("計算完了", "完了通知")

    End Sub

    Private Sub main_sng_pin(ByVal sw As System.IO.StreamWriter, ByVal k As Integer)
        Dim n As Integer = 7
        Dim i As Integer
        Dim j As Integer
        Dim a(n, n + 1) As Double
        Dim ra As Double : Dim rb As Double : Dim rr As Double
        Dim r1 As Double : Dim r2 As Double
        Dim Eg As Double : Dim ng As Double
        Dim Es As Double : Dim ns As Double : Dim alps As Double
        Dim Ec As Double : Dim alpc As Double
        Dim P0 As Double : Dim DT As Double
        Dim C1 As Double : Dim C2 As Double
        Dim xx As Double : Dim xa As Double
        Dim E As Double : Dim nu As Double : Dim alp As Double : Dim TT As Double
        Dim u As Double : Dim sigr As Double : Dim sigt As Double
        Dim ug As Double : Dim sigrg As Double : Dim sigtg As Double
        Dim us2 As Double : Dim sigrs2 As Double : Dim sigts2 As Double
        Dim us1 As Double : Dim sigrs1 As Double : Dim sigts1 As Double
        Dim uc As Double : Dim sigrc As Double : Dim sigtc As Double = 0.0
        Dim dat As String

        ra = aa_inp(k)
        r1 = aa_inp(k) + cc_inp(k)
        r2 = r1 + ta_inp(k)
        rb = bb_inp(k)
        rr = rr_inp(k)
        Eg = Eg_inp(k)
        ng = ng_inp(k)
        Es = Es_inp(k)
        ns = ns_inp(k)
        alps = as_inp(k)
        Ec = Ec_inp(k)
        alpc = ac_inp(k)
        P0 = PP_inp(k) ' Positive : internal pressure 
        DT = TT_inp(k)

        For i = 0 To n
            For j = 0 To n + 1
                a(i, j) = 0.0
        Next j, i

        a(0, 0) = rr : a(0, 1) = 1.0 / rr
        a(1, 0) = Eg / (1.0 + ng) / (1.0 - 2.0 * ng) : a(1, 1) = -Eg / (1.0 + ng) / rb / rb : a(1, 2) = 0.0 : a(1, 3) = -Ec / rb
        a(2, 0) = rb : a(2, 1) = 1.0 / rb : a(2, 2) = -1.0 : a(2, 3) = -Math.Log(rb)
        a(3, 2) = 0.0 : a(3, 3) = Ec / r2 : a(3, 4) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(3, 5) = Es / (1.0 + ns) / r2 / r2
        a(4, 2) = 1.0 : a(4, 3) = Math.Log(r2) : a(4, 4) = -r2 : a(4, 5) = -1.0 / r2
        a(5, 4) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(5, 5) = -Es / (1.0 + ns) / r1 / r1 : a(5, 6) = 0.0 : a(5, 7) = -Ec / r1
        a(6, 4) = r1 : a(6, 5) = 1.0 / r1 : a(6, 6) = -1.0 : a(6, 7) = -Math.Log(r1)
        a(7, 6) = 0.0 : a(7, 7) = Ec / ra

        a(0, 8) = 0.0
        a(1, 8) = 0.0
        a(2, 8) = alpc * DT * (rb - r2)
        a(3, 8) = -Es * alps * DT / (1.0 - ns) * (r2 * r2 - r1 * r1) / 2.0 / r2 / r2
        a(4, 8) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r2 * r2 - r1 * r1) / 2.0 / r2
        a(5, 8) = 0.0
        a(6, 8) = alpc * DT * (r1 - ra)
        a(7, 8) = -P0

        Call MATGJ(n, a)

        C1 = a(0, 8) : C2 = a(1, 8) : xx = rb : xa = rb : E = Eg : nu = ng : alp = 0.0 : TT = 0.0
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        ug = u : sigrg = sigr : sigtg = sigt

        C1 = a(4, 8) : C2 = a(5, 8) : xx = r2 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us2 = u : sigrs2 = sigr : sigts2 = sigt

        C1 = a(4, 8) : C2 = a(5, 8) : xx = r1 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us1 = u : sigrs1 = sigr : sigts1 = sigt

        C1 = a(6, 8) : C2 = a(7, 8) : xx = ra : xa = ra : E = Ec : alp = alpc : TT = DT
        Call CON(C1, C2, xx, xa, E, alp, DT, u, sigr)
        uc = u : sigrc = sigr

        dat = k.ToString()
        dat = dat & "," & IE_inp(k).ToString()
        dat = dat & "," & Eg.ToString("#")
        dat = dat & "," & DT.ToString("#.#")
        dat = dat & "," & sigrc.ToString("#.###")
        dat = dat & "," & sigtc.ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs1 + sigrs2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts1 + sigts2)).ToString("#.###")
        dat = dat & "," & "---"
        dat = dat & "," & "---"
        dat = dat & "," & sigrg.ToString("#.###")
        dat = dat & "," & sigtg.ToString("#.###")
        dat = dat & "," & uc.ToString("#.###")
        dat = dat & "," & ug.ToString("#.###")
        sw.WriteLine(dat)

    End Sub

    Private Sub main_dbl_pin(ByVal sw As System.IO.StreamWriter, ByVal k As Integer)
        Dim n As Integer = 11
        Dim i As Integer
        Dim j As Integer
        Dim a(n, n + 1) As Double
        Dim ra As Double : Dim rb As Double : Dim rr As Double
        Dim r1 As Double : Dim r2 As Double : Dim r3 As Double : Dim r4 As Double
        Dim Eg As Double : Dim ng As Double
        Dim Es As Double : Dim ns As Double : Dim alps As Double
        Dim Ec As Double : Dim alpc As Double
        Dim P0 As Double : Dim DT As Double
        Dim C1 As Double : Dim C2 As Double
        Dim xx As Double : Dim xa As Double
        Dim E As Double : Dim nu As Double : Dim alp As Double : Dim TT As Double
        Dim u As Double : Dim sigr As Double : Dim sigt As Double
        Dim ug As Double : Dim sigrg As Double : Dim sigtg As Double
        Dim us4 As Double : Dim sigrs4 As Double : Dim sigts4 As Double
        Dim us3 As Double : Dim sigrs3 As Double : Dim sigts3 As Double
        Dim us2 As Double : Dim sigrs2 As Double : Dim sigts2 As Double
        Dim us1 As Double : Dim sigrs1 As Double : Dim sigts1 As Double
        Dim uc As Double : Dim sigrc As Double : Dim sigtc As Double = 0.0
        Dim dat As String

        ra = aa_inp(k)
        r1 = ra + cc_inp(k)
        r2 = r1 + ta_inp(k)
        rb = bb_inp(k)
        r4 = rb - cc_inp(k)
        r3 = r4 - tb_inp(k)
        rr = rr_inp(k)
        Eg = Eg_inp(k)
        ng = ng_inp(k)
        Es = Es_inp(k)
        ns = ns_inp(k)
        alps = as_inp(k)
        Ec = Ec_inp(k)
        alpc = ac_inp(k)
        P0 = PP_inp(k) ' Positive : internal pressure 
        DT = TT_inp(k)

        For i = 0 To n
            For j = 0 To n + 1
                a(i, j) = 0.0
        Next j, i

        a(0, 0) = rr : a(0, 1) = 1.0 / rr
        a(1, 0) = Eg / (1.0 + ng) / (1.0 - 2.0 * ng) : a(1, 1) = -Eg / (1.0 + ng) / rb / rb : a(1, 2) = 0.0 : a(1, 3) = -Ec / rb
        a(2, 0) = rb : a(2, 1) = 1.0 / rb : a(2, 2) = -1.0 : a(2, 3) = -Math.Log(rb)
        a(3, 2) = 0.0 : a(3, 3) = Ec / r4 : a(3, 4) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(3, 5) = Es / (1.0 + ns) / r4 / r4
        a(4, 2) = 1.0 : a(4, 3) = Math.Log(r4) : a(4, 4) = -r4 : a(4, 5) = -1.0 / r4
        a(5, 4) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(5, 5) = -Es / (1.0 + ns) / r3 / r3 : a(5, 6) = 0.0 : a(5, 7) = -Ec / r3
        a(6, 4) = r3 : a(6, 5) = 1.0 / r3 : a(6, 6) = -1.0 : a(6, 7) = -Math.Log(r3)
        a(7, 6) = 0.0 : a(7, 7) = Ec / r2 : a(7, 8) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(7, 9) = Es / (1.0 + ns) / r2 / r2
        a(8, 6) = 1.0 : a(8, 7) = Math.Log(r2) : a(8, 8) = -r2 : a(8, 9) = -1.0 / r2
        a(9, 8) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(9, 9) = -Es / (1.0 + ns) / r1 / r1 : a(9, 10) = 0.0 : a(9, 11) = -Ec / r1
        a(10, 8) = r1 : a(10, 9) = 1.0 / r1 : a(10, 10) = -1.0 : a(10, 11) = -Math.Log(r1)
        a(11, 10) = 0.0 : a(11, 11) = Ec / ra

        a(0, 12) = 0.0
        a(1, 12) = 0.0
        a(2, 12) = alpc * DT * (rb - r2)
        a(3, 12) = -Es * alps * DT / (1.0 - ns) * (r4 * r4 - r3 * r3) / 2.0 / r4 / r4
        a(4, 12) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r4 * r4 - r3 * r3) / 2.0 / r4
        a(5, 12) = 0.0
        a(6, 12) = alpc * DT * (r3 - r2)
        a(7, 12) = -Es * alps * DT / (1.0 - ns) * (r2 * r2 - r1 * r1) / 2.0 / r2 / r2
        a(8, 12) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r2 * r2 - r1 * r1) / 2.0 / r2
        a(9, 12) = 0.0
        a(10, 12) = alpc * DT * (r1 - ra)
        a(11, 12) = -P0

        Call MATGJ(n, a)

        C1 = a(0, 12) : C2 = a(1, 12) : xx = rb : xa = rb : E = Eg : nu = ng : alp = 0.0 : TT = 0.0
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        ug = u : sigrg = sigr : sigtg = sigt

        C1 = a(4, 12) : C2 = a(5, 12) : xx = r4 : xa = r3 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us4 = u : sigrs4 = sigr : sigts4 = sigt

        C1 = a(4, 12) : C2 = a(5, 12) : xx = r3 : xa = r3 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us3 = u : sigrs3 = sigr : sigts3 = sigt

        C1 = a(8, 12) : C2 = a(9, 12) : xx = r2 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us2 = u : sigrs2 = sigr : sigts2 = sigt

        C1 = a(8, 12) : C2 = a(9, 12) : xx = r1 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us1 = u : sigrs1 = sigr : sigts1 = sigt

        C1 = a(10, 12) : C2 = a(11, 12) : xx = ra : xa = ra : E = Ec : alp = alpc : TT = DT
        Call CON(C1, C2, xx, xa, E, alp, DT, u, sigr)
        uc = u : sigrc = sigr

        dat = k.ToString()
        dat = dat & "," & IE_inp(k).ToString()
        dat = dat & "," & Eg.ToString("#")
        dat = dat & "," & DT.ToString("#.#")
        dat = dat & "," & sigrc.ToString("#.###")
        dat = dat & "," & sigtc.ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs1 + sigrs2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts1 + sigts2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs3 + sigrs4)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts3 + sigts4)).ToString("#.###")
        dat = dat & "," & sigrg.ToString("#.###")
        dat = dat & "," & sigtg.ToString("#.###")
        dat = dat & "," & uc.ToString("#.###")
        dat = dat & "," & ug.ToString("#.###")
        sw.WriteLine(dat)

    End Sub

    Private Sub main_sng_pex(ByVal sw As System.IO.StreamWriter, ByVal k As Integer)
        Dim n As Integer = 5
        Dim i As Integer
        Dim j As Integer
        Dim a(n, n + 1) As Double
        Dim ra As Double : Dim rb As Double
        Dim r1 As Double : Dim r2 As Double
        Dim Eg As Double = 0.0
        Dim Es As Double : Dim ns As Double : Dim alps As Double
        Dim Ec As Double : Dim nc As Double : Dim alpc As Double
        Dim Pe As Double : Dim DT As Double
        Dim C1 As Double : Dim C2 As Double
        Dim xx As Double : Dim xa As Double
        Dim E As Double : Dim nu As Double : Dim alp As Double : Dim TT As Double
        Dim u As Double : Dim sigr As Double : Dim sigt As Double
        Dim ug As Double : Dim sigrg As Double : Dim sigtg As Double
        Dim us2 As Double : Dim sigrs2 As Double : Dim sigts2 As Double
        Dim us1 As Double : Dim sigrs1 As Double : Dim sigts1 As Double
        Dim uc As Double : Dim sigrc As Double : Dim sigtc As Double = 0.0
        Dim dat As String

        ra = aa_inp(k)
        r1 = aa_inp(k) + cc_inp(k)
        r2 = r1 + ta_inp(k)
        rb = bb_inp(k)
        Es = Es_inp(k)
        ns = ns_inp(k)
        alps = as_inp(k)
        Ec = Ec_inp(k)
        nc = nc_inp(k)
        alpc = ac_inp(k)
        Pe = PP_inp(k) ' Positive : outer pressure 
        DT = TT_inp(k)

        For i = 0 To n
            For j = 0 To n + 1
                a(i, j) = 0.0
        Next j, i

        a(0, 0) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(0, 1) = -Ec / (1.0 + nc) / rb / rb
        a(1, 0) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(1, 1) = -Ec / (1.0 + nc) / r2 / r2 : a(1, 2) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(1, 3) = Es / (1.0 + ns) / r2 / r2
        a(2, 0) = r2 : a(2, 1) = 1.0 / r2 : a(2, 2) = -r2 : a(2, 3) = -1.0 / r2
        a(3, 2) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(3, 3) = -Es / (1.0 + ns) / r1 / r1 : a(3, 4) = -Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(3, 5) = Ec / (1.0 + nc) / r1 / r1
        a(4, 2) = r1 : a(4, 3) = 1.0 / r1 : a(4, 4) = -r1 : a(4, 5) = -1.0 / r1
        a(5, 4) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(5, 5) = -Ec / (1.0 + nc) / ra / ra

        a(0, 6) = -Pe
        a(1, 6) = -Es * alps * DT / (1.0 - ns) * (r2 * r2 - r1 * r1) / 2.0 / r2 / r2
        a(2, 6) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r2 * r2 - r1 * r1) / 2.0 / r2
        a(3, 6) = -Ec * alpc * DT / (1.0 - nc) * (r1 * r1 - ra * ra) / 2.0 / r1 / r1
        a(4, 6) = (1.0 + nc) / (1.0 - nc) * alpc * DT * (r1 * r1 - ra * ra) / 2.0 / r1
        a(5, 6) = 0.0

        Call MATGJ(n, a)

        C1 = a(0, 6) : C2 = a(1, 6) : xx = rb : xa = r2 : E = Ec : nu = nc : alp = alpc : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        ug = u : sigrg = sigr : sigtg = sigt

        C1 = a(2, 6) : C2 = a(3, 6) : xx = r2 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us2 = u : sigrs2 = sigr : sigts2 = sigt

        C1 = a(2, 6) : C2 = a(3, 6) : xx = r1 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us1 = u : sigrs1 = sigr : sigts1 = sigt

        C1 = a(4, 6) : C2 = a(5, 6) : xx = ra : xa = ra : E = Ec : nu = nc : alp = alpc : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        uc = u : sigrc = sigr : sigtc = sigt

        dat = k.ToString()
        dat = dat & "," & IE_inp(k).ToString()
        dat = dat & "," & Eg.ToString("#")
        dat = dat & "," & DT.ToString("#.#")
        dat = dat & "," & sigrc.ToString("#.###")
        dat = dat & "," & sigtc.ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs1 + sigrs2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts1 + sigts2)).ToString("#.###")
        dat = dat & "," & "---"
        dat = dat & "," & "---"
        dat = dat & "," & sigrg.ToString("#.###")
        dat = dat & "," & sigtg.ToString("#.###")
        dat = dat & "," & uc.ToString("#.###")
        dat = dat & "," & ug.ToString("#.###")
        sw.WriteLine(dat)

    End Sub

    Private Sub main_dbl_pex(ByVal sw As System.IO.StreamWriter, ByVal k As Integer)
        Dim n As Integer = 9
        Dim i As Integer
        Dim j As Integer
        Dim a(n, n + 1) As Double
        Dim ra As Double : Dim rb As Double
        Dim r1 As Double : Dim r2 As Double : Dim r3 As Double : Dim r4 As Double
        Dim Eg As Double = 0.0
        Dim Es As Double : Dim ns As Double : Dim alps As Double
        Dim Ec As Double : Dim nc As Double : Dim alpc As Double
        Dim Pe As Double : Dim DT As Double
        Dim C1 As Double : Dim C2 As Double
        Dim xx As Double : Dim xa As Double
        Dim E As Double : Dim nu As Double : Dim alp As Double : Dim TT As Double
        Dim u As Double : Dim sigr As Double : Dim sigt As Double
        Dim ug As Double : Dim sigrg As Double : Dim sigtg As Double
        Dim us4 As Double : Dim sigrs4 As Double : Dim sigts4 As Double
        Dim us3 As Double : Dim sigrs3 As Double : Dim sigts3 As Double
        Dim us2 As Double : Dim sigrs2 As Double : Dim sigts2 As Double
        Dim us1 As Double : Dim sigrs1 As Double : Dim sigts1 As Double
        Dim uc As Double : Dim sigrc As Double : Dim sigtc As Double = 0.0
        Dim dat As String

        ra = aa_inp(k)
        r1 = ra + cc_inp(k)
        r2 = r1 + ta_inp(k)
        rb = bb_inp(k)
        r4 = rb - cc_inp(k)
        r3 = r4 - tb_inp(k)
        Es = Es_inp(k)
        ns = ns_inp(k)
        alps = as_inp(k)
        Ec = Ec_inp(k)
        nc = nc_inp(k)
        alpc = ac_inp(k)
        Pe = PP_inp(k) ' Positive : outer pressure 
        DT = TT_inp(k)

        For i = 0 To n
            For j = 0 To n + 1
                a(i, j) = 0.0
        Next j, i

        a(0, 0) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(0, 1) = -Ec / (1.0 + nc) / rb / rb
        a(1, 0) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(1, 1) = -Ec / (1.0 + nc) / r4 / r4 : a(1, 2) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(1, 3) = Es / (1.0 + ns) / r4 / r4
        a(2, 0) = r4 : a(2, 1) = 1.0 / r4 : a(2, 2) = -r4 : a(2, 3) = -1.0 / r4
        a(3, 2) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(3, 3) = -Es / (1.0 + ns) / r3 / r3 : a(3, 4) = -Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(3, 5) = Ec / (1.0 + nc) / r3 / r3
        a(4, 2) = r3 : a(4, 3) = 1.0 / r3 : a(4, 4) = -r3 : a(4, 5) = -1.0 / r3
        a(5, 4) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(5, 5) = -Ec / (1.0 + nc) / r2 / r2 : a(5, 6) = -Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(5, 7) = Es / (1.0 + ns) / r2 / r2
        a(6, 4) = r2 : a(6, 5) = 1.0 / r2 : a(6, 6) = -r2 : a(6, 7) = -1.0 / r2
        a(7, 6) = Es / (1.0 + ns) / (1.0 - 2.0 * ns) : a(7, 7) = -Es / (1.0 + ns) / r1 / r1 : a(7, 8) = -Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(7, 9) = Ec / (1.0 + nc) / r1 / r1
        a(8, 6) = r1 : a(8, 7) = 1.0 / r1 : a(8, 8) = -r1 : a(8, 9) = -1.0 / r1
        a(9, 8) = Ec / (1.0 + nc) / (1.0 - 2.0 * nc) : a(9, 9) = -Ec / (1.0 + nc) / ra / ra

        a(0, 10) = -Pe
        a(1, 10) = -Es * alps * DT / (1.0 - ns) * (r4 * r4 - r3 * r3) / 2.0 / r4 / r4
        a(2, 10) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r4 * r4 - r3 * r3) / 2.0 / r4
        a(3, 10) = -Ec * alpc * DT / (1.0 - nc) * (r3 * r3 - r2 * r2) / 2.0 / r3 / r3
        a(4, 10) = (1.0 + nc) / (1.0 - nc) * alpc * DT * (r3 * r3 - r2 * r2) / 2.0 / r3
        a(5, 10) = -Es * alps * DT / (1.0 - ns) * (r2 * r2 - r1 * r1) / 2.0 / r2 / r2
        a(6, 10) = (1.0 + ns) / (1.0 - ns) * alps * DT * (r2 * r2 - r1 * r1) / 2.0 / r2
        a(7, 10) = -Ec * alpc * DT / (1.0 - nc) * (r1 * r1 - ra * ra) / 2.0 / r1 / r1
        a(8, 10) = (1.0 + nc) / (1.0 - nc) * alpc * DT * (r1 * r1 - ra * ra) / 2.0 / r1
        a(9, 10) = 0.0

        Call MATGJ(n, a)

        C1 = a(0, 10) : C2 = a(1, 10) : xx = rb : xa = r4 : E = Ec : nu = nc : alp = alpc : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        ug = u : sigrg = sigr : sigtg = sigt

        C1 = a(2, 10) : C2 = a(3, 10) : xx = r4 : xa = r3 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us4 = u : sigrs4 = sigr : sigts4 = sigt

        C1 = a(2, 10) : C2 = a(3, 10) : xx = r3 : xa = r3 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us3 = u : sigrs3 = sigr : sigts3 = sigt

        C1 = a(6, 10) : C2 = a(7, 10) : xx = r2 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us2 = u : sigrs2 = sigr : sigts2 = sigt

        C1 = a(6, 10) : C2 = a(7, 10) : xx = r1 : xa = r1 : E = Es : nu = ns : alp = alps : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        us1 = u : sigrs1 = sigr : sigts1 = sigt

        C1 = a(8, 10) : C2 = a(9, 10) : xx = ra : xa = ra : E = Ec : nu = nc : alp = alpc : TT = DT
        Call ELAS(C1, C2, xx, xa, E, nu, alp, TT, u, sigr, sigt)
        uc = u : sigrc = sigr : sigtc = sigt

        dat = k.ToString()
        dat = dat & "," & IE_inp(k).ToString()
        dat = dat & "," & Eg.ToString("#")
        dat = dat & "," & DT.ToString("#.#")
        dat = dat & "," & sigrc.ToString("#.###")
        dat = dat & "," & sigtc.ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs1 + sigrs2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts1 + sigts2)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigrs3 + sigrs4)).ToString("#.###")
        dat = dat & "," & (0.5 * (sigts3 + sigts4)).ToString("#.###")
        dat = dat & "," & sigrg.ToString("#.###")
        dat = dat & "," & sigtg.ToString("#.###")
        dat = dat & "," & uc.ToString("#.###")
        dat = dat & "," & ug.ToString("#.###")
        sw.WriteLine(dat)

    End Sub

    Private Sub ELAS(ByVal C1 As Double, ByVal C2 As Double, ByVal xx As Double, ByVal xa As Double, ByVal E As Double, ByVal nu As Double, ByVal alp As Double, ByVal TT As Double, _
                     ByRef u As Double, ByRef sigr As Double, ByRef sigt As Double)
        u = (1.0 + nu) / (1.0 - nu) * alp * TT * (xx * xx - xa * xa) / 2.0 / xx + C1 * xx + C2 / xx
        sigr = -E * alp * TT / (1.0 - nu) * (xx * xx - xa * xa) / 2.0 / xx / xx + E / (1.0 + nu) / (1.0 - 2.0 * nu) * C1 - E / (1.0 + nu) * C2 / xx / xx
        sigt = -E * alp * TT / (1.0 - nu) * (xx * xx + xa * xa) / 2.0 / xx / xx + E / (1.0 + nu) / (1.0 - 2.0 * nu) * C1 + E / (1.0 + nu) * C2 / xx / xx
    End Sub

    Private Sub CON(ByVal C1 As Double, ByVal C2 As Double, ByVal xx As Double, ByVal xa As Double, ByVal E As Double, ByVal alp As Double, ByVal TT As Double, ByRef u As Double, ByRef sigr As Double)
        u = C1 + C2 * Math.Log(xx) + alp * TT * (xx - xa)
        sigr = E * C2 / xx
    End Sub

    Private Sub MATGJ(ByVal n As Integer, ByRef a(,) As Double)
        'Gauss-Jordan法による連立一次方程式の解法
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim s As Integer
        Dim p As Double
        Dim d As Double
        Dim max As Double
        Dim dumy As Double

        For k = 0 To n
            max = 0.0
            s = k
            For j = k To n
                If Math.Abs(a(j, k)) > max Then
                    max = Math.Abs(a(j, k))
                    s = j
                End If
            Next j
            For j = 0 To n + 1
                dumy = a(k, j)
                a(k, j) = a(s, j)
                a(s, j) = dumy
            Next j
            p = a(k, k)
            For j = k To n + 1
                a(k, j) = a(k, j) / p
            Next j
            For i = 0 To n
                If i <> k Then
                    d = a(i, k)
                    For j = k To n + 1
                        a(i, j) = a(i, j) - d * a(k, j)
                    Next j
                End If
            Next i
        Next k

    End Sub

End Class