24小时热门版块排行榜

返回列表

本帖产生 1 个程序强帖，点击这里进行查看

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

libralibra

至尊木虫 (著名写手)

骠骑将军

程序强帖: 40
应助: 817 (博后)
金币: 12914.1
红花: 64
帖子: 2238
在线: 287.3小时
虫号: 696514
注册: 2009-02-05
专业: 计算机软件

【答案】应助回帖

★ ★
xzhdty(金币+2): 欢迎常来程序语言讨论 2011-11-11 10:48:28
dubo(金币+1000): 100,明天我测试一下给其他的金币，呵呵，辛苦了！！！ 2011-11-11 23:33:21
ben_ladeng(程序强帖+1): 专家考核存档 2011-11-12 23:28:11

非常抱歉,看错了gx和gx2,帖子无法修改,重贴了

CODE:

 #include 

#include 

#include 

// --------------------------

// load data from xvg file

// --------------------------

void loadData(FILE *fid, float data[797][2])

{

        int i=0;

        while(fscanf(fid,"%f %f\n",&data[i][0],&data[i][1])!=EOF)

                i++;

}

// --------------------------

// write back to xvg file

// --------------------------

void saveData(float data[1501][7])

{

        FILE *fout;

        int i,j;

        fout = fopen("table_CR1_CR1.xvg","w");

        if(fout==NULL)

                return;

        for (i=0;i<1501;++i)

        {

                for (j=0;j<7;++j)

                {

                        fprintf(fout,"%f ",data[i][j]);

                }

                fprintf(fout,"\n");

        }

        fclose(fout);

}

// --------------------------

// zeros vector

// --------------------------

void zeros(float vec[797])

{

        int i;

        for (i=0;i<797;++i)

                vec[i] = 0.0;

}

// --------------------------

// smooth the array with 

// default span 5

// --------------------------

void smooth(float vec[797])

{

        float sum;

        int i,j;

        for (i=0;i<797;++i)

        {

                sum = 0.0;

                if (i<3)

                {

                        for (j=0;j<=i;++j)

                                sum += vec[j];

                        vec[i] = sum/(i+1);

                }

                else if (i>=797-3)

                {

                        for (j=i;j<797;++j)

                                sum += vec[j];

                        vec[i] = sum/(797-i);

                }

                else

                {

                        for (j=i-2;j<=i+2;++j)

                                sum += vec[j];

                        vec[i] = sum/5;

                }

        }

}

// --------------------------

// main entry

// --------------------------

int main(int argc, char *argv[])

{

         FILE *fid;

        float eps = 2.2204e-16;

        float data[797][2], tablenon[1501][7];

        float x[797], dx[797],

                fx[797], fx2[797],

                hx[797], hx2[797],

                gx[797], gx2[797];

        double xfrom = 0.0, xto = 1.592, xtobig = 3.0, jingdu = 0.002, kbt;

        int T = 500, multi = 500, from = 1, to = 797, len = 797;

        int i,j;

         fid = fopen("rdf.xvg","r");

         if(fid==NULL)

                 return -1;

        // read to data 797*2

         loadData(fid,data);

         fclose(fid);

        // make x and dx

        for (i=0;i<797;++i)

        {

                x[i] = i*jingdu;

                dx[i] = data[i][1];

        }

        // update dx

        kbt = T*8.3145/1000;

        smooth(dx);

        for (i=0;i<797;++i)

        {

                dx[i] += eps;

                dx[i] = -1*kbt*log(dx[i]);

        }

        // smooth 10 times

        for        (i=0;i<10;++i)

                smooth(dx);

        zeros(fx);

        zeros(fx2);

        zeros(hx);

        zeros(hx2);

        zeros(gx);

        zeros(gx2);

        for (i=0;i<797;++i)

                gx[i] = dx[i];

        smooth(gx2);

        // first 797 row

        for (i=0;i<797;++i)

        {

                tablenon[i][0] = x[i];

                tablenon[i][1] = fx[i];

                tablenon[i][2] = fx2[i];

                tablenon[i][3] = gx[i];

                tablenon[i][4] = gx2[i];

                tablenon[i][5] = hx[i];

                tablenon[i][6] = hx2[i];

        }

        // the rest

        for (i=797;i<1501;++i)

        {

                tablenon[i][0] = xto+(i+1-797)*jingdu;

                for (j=1;j<7;++j)

                        tablenon[i][j] = 0;

        }

        saveData(tablenon);

        system("PAUSE");

        return 0;

}