C++仿真微带电路

//微带线结构，激励为平面波激励。采用shape函数限制平面波在微带线和地之间
//高斯脉冲
//介质为微带结构
#define _CRT_SECURE_NO_DEPRECATE
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define IE 62
#define JE 120
#define KE 14

#define ia 14
#define ja 15
#define ka 7
#define NFREQS 3
#define ktop 2
void main()
{
    double dx[IE][JE][KE], dy[IE][JE][KE], dz[IE][JE][KE];
    double ex[IE][JE][KE], ey[IE][JE][KE], ez[IE][JE][KE];
    double hx[IE][JE][KE], hy[IE][JE][KE], hz[IE][JE][KE];
    //double ix[IE][JE][KE], iy[IE][JE][KE], iz[IE][JE][KE];
    double gax[IE][JE][KE], gay[IE][JE][KE], gaz[IE][JE][KE];
    double gbx[IE][JE][KE], gby[IE][JE][KE], gbz[IE][JE][KE];
    int   n, i, j, k, ic, jc, kc, nsteps, npml;
    double ddz, ra_x, ra_y, dt, T, epsz, muz, pi;
    int ib, jb, kb;
    double  xn, xxn, xnum, xd;
    double t0, spread, pulse;
    double ez_inc[JE], hx_inc[JE];
    double ez_inc_low_m1, ez_inc_low_m2, ez_inc_high_m1, ez_inc_high_m2;

    double idxl[ia][JE][KE], idxh[ia][JE][KE];
    double ihxl[ia][JE][KE], ihxh[ia][JE][KE];
    double idyl[IE][ja][KE], idyh[IE][ja][KE];
    double ihyl[IE][ja][KE], ihyh[IE][ja][KE];
    double idzl[IE][JE][ka], idzh[IE][JE][ka];
    double ihzl[IE][JE][ka], ihzh[IE][JE][ka];
    int ixh, jyh, kzh;

    double gi1[IE], gi2[IE], gi3[IE];
    double gj1[JE], gj2[JE], gj3[JE];
    double gk1[KE], gk2[KE], gk3[KE];
    double fi1[IE], fi2[IE], fi3[IE];
    double fj1[JE], fj2[JE], fj3[JE];
    double fk1[KE], fk2[KE], fk3[KE];

    double curl_h, curl_e;

//    int ii, jj, kk;
//    double dist, xdist, ydist, zdist;

    int istart, iend, k_ref, i_ref;
    int j_patch_st, j_patch_end, j_ref;
    int i_patch_st, i_patch_end;
    double eps_sub;
    double shape[IE][KE], half_wv;
    FILE *fp;

    ic = IE / 2;
    jc = JE / 2;
    kc = KE / 2;
    //PML边界空间大小
    ib = IE - ia - 1;
    jb = JE - ja - 1;
    kb = KE - ka - 1;

    i_patch_st = ia + 5;
    i_patch_end = i_patch_st + 31;
    j_patch_end = jb - 5;
    j_patch_st = j_patch_end - 39;
    j_ref = j_patch_st - 30;
    k_ref = ktop - 1;

    pi = 3.14159;
    npml = 7;
    epsz = 8.8e-12;         //真空中的介电常数
    muz = 4 * pi*1.0e-7;    //真空中的磁导率
    ddz = 0.265e-3;             // 空间步长，单个波长中取20个网格
    dt = ddz / 6e8;           //时间步长
    ra_x = 0.6812;
    ra_y = 0.6625;            //x和y方向的网格因子
    eps_sub = 2.2;

    ez_inc_low_m1 = 0.0;
    ez_inc_low_m2 = 0.0;
    ez_inc_high_m1 = 0.0;
    ez_inc_high_m2 = 0.0;

    // Initialize the arrays
    for (j = 0; j < JE; j++){
        ez_inc[j] = 0.0;
        hx_inc[j] = 0.0;
        for (k = 0; k < KE; k++){
            for (i = 0; i < IE; i++){
                ex[i][j][k] = 0.0;
                ey[i][j][k] = 0.0;
                ez[i][j][k] = 0.0;
                dx[i][j][k] = 0.0;
                dy[i][j][k] = 0.0;
                dz[i][j][k] = 0.0;
                hx[i][j][k] = 0.0;
                hy[i][j][k] = 0.0;
                hz[i][j][k] = 0.0;
                gax[i][j][k] = 1.0;
                gay[i][j][k] = 1.0;
                gaz[i][j][k] = 1.0;
                gbx[i][j][k] = 0.0;
                gby[i][j][k] = 0.0;
                gbz[i][j][k] = 0.0;
            }
        }
    }
    //定义基板
    for (j = 0; j < JE; j++)
    {
        for (i = 0; i < IE; i++)
        {
            for (k = 0; k <= ktop; k++)
            {
                gax[i][j][k] = 1.0 / eps_sub;
                gay[i][j][k] = 1.0 / eps_sub;
                gaz[i][j][k] = 1.0 / eps_sub;
            }
        }
    }
    //底部加一层金属片
    for (i = 1; i < IE - 1; i++)
    {
        for (j = 1; j < JE - 1; j++)
        {
            k = 0;
            gax[i][j][k] = 0.0;
            gay[i][j][k] = 0.0;
        }
    }
    istart = ia + 13;
    iend = istart + 6;
    i_ref = istart + (iend - istart) / 2;
    half_wv = (iend - istart) / 2;

    for (i = istart; i <= iend; i++)
    {
        for (k = 0; k <= ktop; k++)
        {
            shape[i][k] = 1.0;
        }
    }

    //在K=ktop+1上增加导体导线
    for (j = 1; j <= j_patch_st; j++)
    {
        for (i = istart; i < iend - 1; i++)
        {
            gax[i][j][ktop + 1] = 0.0;
        }
        for (i = istart; i < iend - 1; i++)
        {
            gay[i][j][ktop + 1] = 0.0;
        }
    }
    //在ktop上加矩形贴片
    for (j = j_patch_st; j <= j_patch_end; j++)
    {
        for (i = ia + 1; i <= ib - 1; i++)
        {
            gax[i][j][ktop + 1] = 0.0;
        }
        for (i = ia + 1; i <= ib - 1; i++)
        {
            gay[i][j][ktop + 1] = 0.0;
        }
    }
    for (i = 0; i < ia; i++)
    {
        for (j = 0; j < JE; j++)
        {
            for (k = 0; k < KE; k++)
            {
                idxl[i][j][k] = 0.0;
                idxh[i][j][k] = 0.0;
                ihxl[i][j][k] = 0.0;
                ihxh[i][j][k] = 0.0;
            }
        }
    }
    for (i = 0; i < IE; i++)
    {
        for (j = 0; j < ja; j++)
        {
            for (k = 0; k < KE; k++)
            {
                idyl[i][j][k] = 0.0;
                idyh[i][j][k] = 0.0;
                ihyl[i][j][k] = 0.0;
                ihyh[i][j][k] = 0.0;
            }
        }
    }
    for (i = 0; i < IE; i++)
    {
        for (j = 0; j < JE; j++)
        {
            for (k = 0; k < ka; k++)
            {
                idzl[i][j][k] = 0.0;
                idzh[i][j][k] = 0.0;
                ihzl[i][j][k] = 0.0;
                ihzh[i][j][k] = 0.0;
            }
        }
    }

    //Boundry Conditions

    for (i = 0; i < IE; i++)
    {
        gi1[i] = 0.0;
        fi1[i] = 0.0;
        gi2[i] = 1.0;
        fi2[i] = 1.0;
        gi3[i] = 1.0;
        fi3[i] = 1.0;
    }
    for (j = 0; j < JE; j++)
    {
        gj1[j] = 0.0;
        fj1[j] = 0.0;
        gj2[j] = 1.0;
        fj2[j] = 1.0;
        gj3[j] = 1.0;
        fj3[j] = 1.0;
    }
    for (k= 0; k < KE; k++)
    {
        gk1[k] = 0.0;
        fk1[k] = 0.0;
        gk2[k] = 1.0;
        fk2[k] = 1.0;
        gk3[k] = 1.0;
        fk3[k] = 1.0;
    }
    //PML边界设置
    for (i = 0; i < npml; i++)
    {
        xnum = npml - i;
        xd = npml;
        xxn = xnum / xd;
        xn = 0.33*pow(xxn, 3.0);
        fi1[i] = xn;
        fi1[IE - 1 - i] = xn;
        gi2[i] = 1.0 / (1.0 + xn);
        gi2[IE - i - 1] = 1.0 / (1.0 + xn);
        gi3[i] = (1.0 - xn) / (1.0 + xn);
        gi3[IE - i - 1] = (1.0 - xn) / (1.0 + xn);

        xxn = (xnum - 0.5) / xd;
        xn = 0.25*pow(xxn, 3.0);
        gi1[i] = xn;
        gi1[IE - i - 2] = xn;
        fi2[i] = 1.0 / (1.0 + xn);
        fi2[IE - 2 - i] = 1.0 / (1.0 + xn);
        fi3[i] = (1.0 - xn) / (1.0 + xn);
        fi3[IE - i - 2] = (1.0 - xn) / (1.0 + xn);
    }
    
    for (j = 0; j < npml; j++)
    {
        xnum = npml - j;
        xd = npml;
        xxn = xnum / xd;
        xn = 0.33*pow(xxn, 3.0);
        fj1[j] = xn;
        fj1[JE - 1 - j] = xn;
        gj2[j] = 1.0 / (1.0 + xn);
        gj2[JE - j - 1] = 1.0 / (1.0 + xn);
        gj3[j] = (1.0 - xn) / (1.0 + xn);
        gj3[JE - j - 1] = (1.0 - xn) / (1.0 + xn);

        xxn = (xnum - 0.5) / xd;
        xn = 0.33*pow(xxn, 3.0);
        gj1[j] = xn;
        gj1[JE - j - 2] = xn;
        fj2[j] = 1.0 / (1.0 + xn);
        fj2[JE - 2 - j] = 1.0 / (1.0 + xn);
        fj3[j] = (1.0 - xn) / (1.0 + xn);
        fj3[JE - 2 - j] = (1.0 - xn) / (1.0 + xn);
    }
    //在K方向上，只有向上的方向设置了PML
    //向下的方向由于金属板的存在而不需要另设边界
    for (k = 0; k < npml; k++)
    {
        xnum = npml - k;
        xd = npml;
        xxn = xnum / xd;
        xn = 0.33*pow(xxn, 3.0);
        //fj1[j] = xn;
        fk1[KE - 1 - k] = xn;
        //gj2[j] = 1.0 / (1.0 + xn);
        gk2[KE - k - 1] = 1.0 / (1.0 + xn);
        //gj3[j] = (1.0 - xn) / (1.0 + xn);
        gk3[KE - k - 1] = (1.0 - xn) / (1.0 + xn);

        xxn = (xnum - 0.5) / xd;
        xn = 0.33*pow(xxn, 3.0);
        //gj1[j] = xn;
        gk1[KE - k - 2] = xn;
        //fj2[j] = 1.0 / (1.0 + xn);
        fk2[KE - 2 - k] = 1.0 / (1.0 + xn);
        //fj3[j] = (1.0 - xn) / (1.0 + xn);
        fk3[KE - 2 - k] = (1.0 - xn) / (1.0 + xn);
    }

    t0 = 150.0;
    spread = 25;
    T = 0;
    nsteps = 1;

    while (nsteps > 0)
    {
        printf("nsteps -->");
        scanf("%d", &nsteps);
        printf("%d\n", nsteps);

        for (n = 1; n < nsteps; n++)
        {
            T = T + 1;

            //Start of the Main FDTD loop

            //一维高斯脉冲激励用于平面波拓展
            for (j = 1; j < JE; j++)
            {
                ez_inc[j] = gj3[j] * ez_inc[j] + gj2[j] * (0.5*ra_y / eps_sub)*(hx_inc[j - 1] - hx_inc[j]);
            }
            //Source
            pulse = exp(-0.5*pow((t0 - T) / spread, 2.0));
            ez_inc[ja - 2] = pulse;

            //Boundary condition for the incident buffer
            //一维吸收边界
            ez_inc[0] = ez_inc_low_m2;
            ez_inc_low_m2 = ez_inc_low_m1;
            ez_inc_low_m1 = ez_inc[1];
            ez_inc[JE - 1] = ez_inc_high_m2;
            ez_inc_high_m2 = ez_inc_high_m1;
            ez_inc_high_m1 = ez_inc[JE - 2];

            //Calculate the Dx field
            for (i = 1; i < ia; i++)
            {
                for (j = 1; j < JE; j++)
                {
                    for (k = 1; k < KE; k++)
                    {
                        curl_h = (ra_y* (hz[i][j][k] - hz[i][j - 1][k]) - hy[i][j][k] + hy[i][j][k - 1]);
                        idxl[i][j][k] = idxl[i][j][k] + curl_h;
                        dx[i][j][k] = gj3[j] * gk3[k] * dx[i][j][k] + gj2[j] * gk2[k] * 0.5*(curl_h + gi1[i] * idxl[i][j][k]);
                    }
                }
            }

            for (i = ia; i <= ib; i++)
            {
                for (j = 1; j < JE; j++)
                {
                    for (k = 1; k < KE; k++)
                    {
                        curl_h = (ra_y*(hz[i][j][k] - hz[i][j - 1][k]) - hy[i][j][k] + hy[i][j][k - 1]);
                        dx[i][j][k] = gj3[j] * gk3[k] * dx[i][j][k] + gj2[j] * gk2[k] * 0.5*curl_h;
                    }
                }
            }

            for (i = ib+1; i < IE; i++)
            {
                ixh = i - ib - 1;
                for (j = 1; j < JE; j++)
                {
                    for (k = 1; k < KE; k++)
                    {
                        curl_h = (ra_y* (hz[i][j][k] - hz[i][j - 1][k]) - hy[i][j][k] + hy[i][j][k - 1]);
                        idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
                        dx[i][j][k] = gj3[j] * gk3[k] * dx[i][j][k] + gj2[j] * gk2[k] * 0.5*(curl_h + gi1[i] * idxh[i][j][k]);
                    }
                }
            }

            //Calculate the Dy field

            for (i = 1; i < IE; i++)
            {
                for (j = 1; j < ja; j++)
                {
                    for (k = 1; k < KE; k++)
                    {
                        curl_h = (-ra_x* (hz[i][j][k] - hz[i-1][j ][k]) + hx[i][j][k] - hx[i][j][k - 1]);
                        idyl[i][j][k] = idyl[i][j][k] + curl_h;
                        dy[i][j][k] = gi3[i] * gk3[k] * dy[i][j][k] + gi2[i] * gk2[k] * 0.5*(curl_h + gj1[j] * idyl[i][j][k]);
                    }
                }
            }

            for (i =1; i <IE; i++)
            {
                for (j = ja; j <= jb; j++)
                {
                    for (k = 1; k < KE; k++)
                    {
                        curl_h = (-ra_x*(hz[i][j][k] - hz[i-1][j][k]) + hx[i][j][k] - hx[i][j][k - 1]);
                        dy[i][j][k] = gi3[i] * gk3[k] * dy[i][j][k] + gi2[i] * gk2[k] * 0.5*curl_h;
                    }
                }
            }
            for (i = 1; i < IE; i++){
                for (j = jb + 1; j < JE; j++){
                    jyh = j - jb - 1;
                    for (k = 1; k < KE; k++){
                        curl_h = (hx[i][j][k] - hx[i][j][k - 1] - ra_x*(hz[i][j][k] + hz[i - 1][j][k]));
                        idyh[i][jyh][k] = idyh[i][jyh][k] + curl_h;
                        dy[i][j][k] = gi3[i] * gk3[k] * dy[i][j][k] + gi2[i] * gk2[k] * 0.5*(curl_h + gj1[j] * idyh[i][jyh][k]);
                    }
                }
            }
            
            //Calculate the Dz field

            for (i = 1; i < IE; i++)
            {
                for (j = 1; j < JE; j++)
                {
                    for (k = 0; k < ka; k++)
                    {
                        curl_h = (ra_x*(hy[i][j][k] - hy[i - 1][j][k]) - ra_y*(hx[i][j][k] + hx[i][j - 1][k]));
                        idzl[i][j][k] = idzl[i][j][k] + curl_h;
                        dz[i][j][k] = gi3[i] * gj3[j] * dz[i][j][k] + gi2[i] * gj2[j] * 0.5*(curl_h + gk1[k] * idzl[i][j][k]);
                    }
                }
            }

            for (i = 1; i < IE; i++)
            {
                for (j = 1; j < JE; j++)
                {
                    for (k = ka; k <= kb; k++)
                    {
                        curl_h = (ra_x*(hy[i][j][k] - hy[i - 1][j][k]) - ra_y*(hx[i][j][k] + hx[i][j - 1][k]));
                        dz[i][j][k] = gi3[i] * gj3[j] * dz[i][j][k] + gi2[i] * gj2[j] * 0.5*curl_h;
                    }
                }
            }

            for (i = 1; i < IE; i++)
            {
                for (j = 1; j < JE; j++)
                {
                    for (k = kb + 1; k < KE; k++)
                    {
                        kzh = k - kb - 1;
                        curl_h = (ra_x*(hy[i][j][k] - hy[i - 1][j][k]) - ra_y*(hx[i][j][k] + hx[i][j - 1][k]));
                        idzh[i][j][kzh] = idzh[i][j][kzh] + curl_h;
                        dz[i][j][k] = gi3[i] * gj3[j] * dz[i][j][k] + gi2[i] * gj2[j] * 0.5*(curl_h + gk1[k] * idzh[i][j][kzh]);
                    }
                }
            }

            //修正 Dz
            for (i = istart; i <= iend; i++)
            {
                for (k = 0; k <= ktop; k++)
                {
                    dz[i][ja][k] = dz[i][ja][k] + (0.5 / eps_sub)*shape[i][k] * hx_inc[ja - 1];
                }

                // Calculate the E from D field
                for (i = 1; i < IE - 1; i++)
                {
                    for (j = 1; j < JE - 1; j++)
                    {
                        for (k = 1; k < KE - 1; k++)
                        {
                            ex[i][j][k] = gax[i][j][k] * (dx[i][j][k]);
                            ey[i][j][k] = gay[i][j][k] * (dy[i][j][k]);
                            ez[i][j][k] = gaz[i][j][k] * (dz[i][j][k]);
                        }
                    }
                }

                // Calculate the incident field
                //一维高斯脉冲激励用于平面波扩展，计算磁场
                for (j = 0; j < JE - 1; j++)
                {
                    hx_inc[j] = fj3[j] * hx_inc[j] + 0.5*fj2[j] * (ez_inc[j] - ez_inc[j + 1]);
                }

                //Calculate the Hx field

                for (i = 0; i < ia; i++)
                {
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ey[i][j][k + 1] - ey[i][j][k] - ra_y*(ez[i][j + 1][k] + ez[i][j][k]));
                            ihxl[i][j][k] = ihxl[i][j][k] + curl_e;
                            hx[i][j][k] = fj3[j] * fk3[k] * hx[i][j][k] + fj2[j] * fk2[k] * 0.5*(curl_e + fi1[i] * ihxl[i][j][k]);
                        }
                    }
                }

                for (i = ia; i <= ib; i++)
                {
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ey[i][j][k + 1] - ey[i][j][k] - ra_y*(ez[i][j + 1][k] + ez[i][j][k]));
                            hx[i][j][k] = fj3[j] * fk3[k] * hx[i][j][k] + fj2[j] * fk2[k] * 0.5*curl_e;
                        }
                    }
                }

                for (i = ib + 1; i < IE; i++)
                {
                    ixh = i - ib - 1;
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ey[i][j][k + 1] - ey[i][j][k] - ra_y*(ez[i][j + 1][k] + ez[i][j][k]));
                            ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
                            hx[i][j][k] = fj3[j] * fk3[k] * hx[i][j][k] + fj2[j] * fk2[k] * 0.5*(curl_e + fi1[i] * ihxh[ixh][j][k]);
                        }
                    }
                }

                // 修正 Hx
                for (i = istart; i <= iend; i++)
                {
                    for (k = 0; k <= ktop; k++)
                    {
                        hx[i][ja - 1][k] = hx[i][ja - 1][k] + (0.5 / eps_sub)*shape[i][k] * ez_inc[ja];
                    }
                }

                // Calculate the Hy field

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = 0; j < ja; j++)
                    {
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ra_x*(ez[i + 1][j][k] - ez[i][j][k]) - ex[i][j][k + 1] + ex[i][j][k]);
                            ihyl[i][j][k] = ihyl[i][j][k] + curl_e;
                            hy[i][j][k] = fi3[i] * fk3[k] * hy[i][j][k] + fi2[i] * fk3[k] * 0.5*(curl_e + fj1[j] * ihyl[i][j][k]);
                        }
                    }
                }

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = ja; j <= jb; j++)
                    {
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ra_x*(ez[i + 1][j][k] - ez[i][j][k]) - ex[i][j][k + 1] + ex[i][j][k]);
                            hy[i][j][k] = fi3[i] * fk3[k] * hy[i][j][k] + fi2[i] * fk3[k] * 0.5*curl_e;
                        }
                    }
                }

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = jb + 1; j < JE; j++)
                    {
                        jyh = j - jb - 1;
                        for (k = 0; k < KE - 1; k++)
                        {
                            curl_e = (ra_x*(ez[i + 1][j][k] - ez[i][j][k]) - ex[i][j][k + 1] + ex[i][j][k]);
                            ihyh[i][jyh][k] = ihyh[i][jyh][k] + curl_e;
                            hy[i][j][k] = fi3[i] * fk3[k] * hy[i][j][k] + fi2[i] * fk3[k] * 0.5*(curl_e + fj1[j] * ihyh[i][jyh][k]);
                        }
                    }
                }

                // Calculate the Hz field

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = 0; k < ka; k++)
                        {
                            curl_e = (ra_y*(ex[i][j + 1][k] - ex[i][j][k]) - ra_x*(ey[i + 1][j][k] + ey[i][j][k]));
                            ihzl[i][j][k] = ihzl[i][j][k] + curl_e;
                            hz[i][j][k] = fi3[i] * fj3[j] * hz[i][j][k] + fi2[i] * fj2[j] * 0.5*(curl_e + fk1[k] * ihzl[i][j][k]);
                        }
                    }
                }

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = ka; k <= kb; k++)
                        {
                            curl_e = (ra_y*(ex[i][j + 1][k] - ex[i][j][k]) - ra_x*(ey[i + 1][j][k] + ey[i][j][k]));
                            hz[i][j][k] = fi3[i] * fj3[j] * hz[i][j][k] + fi2[i] * fj2[j] * 0.5*curl_e;
                        }
                    }
                }

                for (i = 0; i < IE - 1; i++)
                {
                    for (j = 0; j < JE - 1; j++)
                    {
                        for (k = kb + 1; k < KE; k++)
                        {
                            kzh = k - kb - 1;
                            curl_e = (ra_y*(ex[i][j + 1][k] - ex[i][j][k]) - ra_x*(ey[i + 1][j][k] + ey[i][j][k]));
                            ihzh[i][j][kzh] = ihzh[i][j][kzh] + curl_e;
                            hz[i][j][k] = fi3[i] * fj3[j] * hz[i][j][k] + fi2[i] * fj2[j] * 0.5*(curl_e + fk1[k] * ihzh[i][j][kzh]);
                        }
                    }
                }
            }

            //End of the main FDTD loop
            //Write the E field out to a file "Ez"
            fp = fopen("Ez.txt", "w");
            for (j = 0; j < JE; j++)
            {
                for (i = 0; i < IE; i++)
                {
                    fprintf(fp, "%9.6f", ez[i][j][ktop]);
                }
                fprintf(fp, "\n");
            }
            fclose(fp);

            printf("T=%4.0f \n",T);
        }
    }
}