【C言語】古典的超伝導体のギャップ方程式からTcを求めてみた

概要

　古典的超伝導体のギャップ方程式からTcを求めた。

動作環境

• Windows10(64bit)
• Python 3.7.7

アルゴリズム

アルゴリズムは二分法を用いた。

プログラム

/* bcs.c */
#include <stdio.h>
#include <math.h>

/* ↓↓↓↓nを大きくすると精度が上がる*/
const int n = 100;
double y;

/*積分する関数 f(x,y)*/
double f(double k_x, double k_y, double kb_T,double myu)
{
double g = 2.5;
double t = 0.1;
double a = 1.0;
double epsilon =-2.0*t*(cos(k_x*a)+cos(k_y*a))-myu;
double E =sqrt(epsilon*epsilon);
return (g/(2.0*E))*tanh(E/(2.0*kb_T));
}

/* 積分の範囲 x0からx_max  yy0からy_max*/
const double x0 = -M_PI;
const double x_max = M_PI;
const double yy0 = -M_PI;
const double y_max = M_PI;

/* x方向の積分 */
double simpson_1(double y, double kb_T,double myu){
    double delta = (x_max-x0)/n;
    double a;
    a = f(x0,y,kb_T,myu);
    a += f(x0+delta*n,y,kb_T,myu);
    int i;
    for(i=1;i<n;i+=2){
        a += 4.0*f(x0+delta*i,y,kb_T,myu);
    }
    for(i=2;i<n;i+=2){
        a += 2.0*f(x0+delta*i,y,kb_T,myu);
    }
    return a*delta/3.0;
}

/* y方向の積分 */
double simpson_2(double half,double myu){
    double delta=(y_max-yy0)/n;
    double a;
    double kb_T = half;
    a = simpson_1(yy0,kb_T,myu);
    a += simpson_1(yy0+delta*n,kb_T,myu);
    int i;
    for(i=1;i<n;i+=2){
        a += 4.0*simpson_1(yy0+delta*i,kb_T,myu);
    }
    for(i=2;i<n;i+=2){
        a += 2.0*simpson_1(yy0+delta*i,kb_T,myu);
    }
    
    return a*delta/3.0;
}

int main(){
int max=1000;
double myu =-0.3;
int dmyu =600;
double array_Tc[dmyu];
double array_myu[dmyu];
double left, half, right;
double residue;

for(int j=0; j<dmyu;j++) #二分法
{
    left=0.0;
    right=20.0;
    for(int k=0; k<max; k++)
    {
        half=(left+right)/2;
        residue=1.0-simpson_2(half,myu);
        if(residue>0)
        {
        right=half;
        }else{
        left=half;
        }
    }
    array_myu[j] = myu; 
    array_Tc[j] = half;
    printf("%d",j);
    myu += 0.001;
}
FILE *fp; #計算結果をtxtファイルに書き出す
fp = fopen("data_c20.txt","w");
for (int l = 0; l <= dmyu; l++)
{
fprintf(fp,"%f %.20f\n",array_myu[l],array_Tc[l]);
}
fclose(fp);
return 0;
}