我需要igrf的C语言代码

以下是一个基于IGRF模型的C语言代码：

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846

int main()
{
    double r, theta, phi; //输入的位置信息
    double B_r, B_theta, B_phi; //计算得到的磁场信息
    double a[13][13], b[13][13]; //IGRF模型系数
    int n, m; //IGRF模型阶数

    //读取位置信息
    printf("Please input the location information:\n");
    printf("r (km): ");
    scanf("%lf", &r);
    printf("theta (deg): ");
    scanf("%lf", &theta);
    printf("phi (deg): ");
    scanf("%lf", &phi);
    theta = theta * PI / 180.0;
    phi = phi * PI / 180.0;

    //读取IGRF模型系数
    printf("Please input the IGRF model order:\n");
    printf("n: ");
    scanf("%d", &n);
    printf("m: ");
    scanf("%d", &m);
    printf("Please input the IGRF model coefficients:\n");
    for(int i=0; i<=n; i++)
    {
        for(int j=0; j<=m; j++)
        {
            scanf("%lf %lf", &a[i][j], &b[i][j]);
        }
    }

    //计算磁场信息
    double x = r * cos(theta);
    double y = r * sin(theta) * cos(phi);
    double z = r * sin(theta) * sin(phi);
    B_r = 0;
    B_theta = 0;
    B_phi = 0;
    for(int i=1; i<=n; i++)
    {
        for(int j=0; j<=i; j++)
        {
            double P = legendre(i, j, sin(theta));
            double cos_m_phi = cos(j * phi);
            double sin_m_phi = sin(j * phi);
            B_r += (a[i][j] * cos_m_phi + b[i][j] * sin_m_phi) * pow(r, -i-1) * (i+1) * P;
            B_theta -= (a[i][j] * cos_m_phi + b[i][j] * sin_m_phi) * pow(r, -i) * P;
            B_phi -= j * (a[i][j] * sin_m_phi - b[i][j] * cos_m_phi) * pow(r, -i) * P;
        }
    }

    //输出磁场信息
    printf("The magnetic field information at (%.2f, %.2f, %.2f) is:\n", r, theta, phi);
    printf("B_r = %.6f (nT)\n", B_r);
    printf("B_theta = %.6f (nT)\n", B_theta);
    printf("B_phi = %.6f (nT)\n", B_phi);

    return 0;
}

double legendre(int n, int m, double x)
{
    double pmm = 1.0;
    if(m > 0)
    {
        double somx2 = sqrt((1.0 - x) * (1.0 + x));
        double fact = 1.0;
        for(int i=1; i<=m; i++)
        {
            pmm *= -fact * somx2;
            fact += 2.0;
        }
    }
    if(n == m)
    {
        return pmm;
    }
    double pmmp1 = x * (2 * m + 1) * pmm;
    if(n == m + 1)
    {
        return pmmp1;
    }
    double pll = 0.0;
    for(int l=m+2; l<=n; l++)
    {
        pll = ((2.0 * l - 1.0) * x * pmmp1 - (l + m - 1.0) * pmm) / (l - m);
        pmm = pmmp1;
        pmmp1 = pll;
    }
    return pll;
}

该代码根据用户输入的位置信息和IGRF模型系数计算出该位置处的地球磁场信息，包括径向、极向和方位向三个分量。其中，`legendre`函数用于计算勒让德多项式。用户需要根据实际需求修改输入输出格式、IGRF模型系数和计算方法。