gnuplotやその他多くのプログラムでは標準関数に近い関数だけど，
これをGSLを利用してやってみようと思いました．以下実装例です．
ループごとにファイル保存なんてことをしているので，実行速度は非常に低速です．
CPUには優しいプログラム（ワライ）です．
今後のプログラムは配列・構造体を導入したりとかして，より速いものにします．
#include <stdio.h>
#include <math.h> //BSD系のC標準数学ライブラリ．
#include <gsl/gsl_sf.h> //特殊関数のヘッダファイル
#include <gsl/gsl_math.h> //数学の基本的な関数のヘッダファイル
#include <gsl/gsl_complex.h> //複素数に関する宣言や定義のヘッダファイル
#include <gsl/gsl_complex_math.h>//複素数を扱った関数のヘッダファイル
int main(void)
{
double theta = 0; //角度theta．初期値を一応設定．
double phi =0; //角度phi．初期値を一応設定．
double YY; //球面調和関数Y．spherical harmonics
//球面調和関数の直交座標系変換した値．
double xx_sh;
double yy_sh;
double zz_sh;
int l=3;
int m=0; //lとm．量子数に対応ですね．
//途中計算で作った変数．数式はTeX記法でやってます．わからない人は適当に調べてね！
//球面調和関数は，Wikipediaによれば次の式．
//Y_{l}^{m}(\theta, \phi) = (-1)^{(m+|m|)/2}
//\sqrt{ \frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!} \,} P_l^{|m|}(\cos\theta) \, e^{im\phi}
//これを分解して見通しをよくする．
//(-1)^{(m+|m|)/2}に対応
double minus1_mm_2;
//\frac{2l+1}{4\pi}に対応
double frac_2l1_4pi;
//\frac{(l-|m|)!}{(l+|m|)!}における，分子に対応
double l_minus_abs_m;
//\frac{(l-|m|)!}{(l+|m|)!}における，分母に対応
double l_plus_abs_m;
//(\cos\theta )\exp (im\phi)に対応
double cos_theta;
double exp_imphi;
int i,j; //ループ用変数．
int stop_i=100; //ループの停止値．まあ設定する必要もない気がするけど作った…．
FILE *outputdata; //データ保存用ストリーム．
FILE *outputdata_xyz; //データ保存用ストリーム．
outputdata = fopen("outputdata.txt", "w"); //データ保存ファイル，オープン
fprintf(outputdata, "#theta\tphi\tYY\n");
outputdata_xyz = fopen("outputdata_xyz.txt", "w"); //データ保存ファイル，オープン
fprintf(outputdata_xyz, "#xx_sh\tyy_sh\tzz_sh\n");
outputdata = fopen("outputdata.txt", "a");
outputdata_xyz = fopen("outputdata_xyz.txt", "a"); //データ保存ファイル，オープン
//計算用forループ．2重になっている．
for (i=1 ; i<=stop_i ; i++){
phi = phi + 0.0628;
for (j=1 ; j<=stop_i ; j++){
theta = theta + 0.0628;
minus1_mm_2 = pow(-1, (m+abs(m))/2);
frac_2l1_4pi = (2*l +1)/(M_PI_4);
l_minus_abs_m = gsl_sf_fact(l-abs(m));
l_plus_abs_m = gsl_sf_fact(l+abs(m));
//cos_theta_exp_imphi = cos(theta)*gsl_expm1();
cos_theta = cos(theta);
exp_imphi = GSL_REAL(gsl_complex_exp(gsl_complex_rect(0, m*phi)))+GSL_IMAG(gsl_complex_exp(gsl_complex_rect(0, m*phi)));
//複素数の扱いで転けた（2011年2月15日，4:49）．寝ます．
//本番の計算．球面調和関数の計算．
YY=frac_2l1_4pi*sqrt(frac_2l1_4pi*(l_minus_abs_m/l_plus_abs_m))*gsl_sf_legendre_Plm(l,m,cos_theta)*exp_imphi;
//↓極座標から直交座標系に変換する．
xx_sh = YY * sin(theta) * cos(phi);
yy_sh = YY * sin(theta) * sin(phi);
zz_sh = YY * cos(theta);
//データ保存用の部分↓
fprintf(outputdata, "%f\t%f\t%f\n", theta,phi,YY);
fprintf(outputdata_xyz, "%f\t%f\t%f\n", xx_sh,yy_sh,zz_sh);
printf("iが%d回目，jが%d回目の計算結果：Y(%f,%f) = %.18e\n", i, j, theta, phi, YY);
}
theta = 0;//ここで初期化しなければthetaの値がめちゃめちゃに大きくなる．
}
fclose(outputdata); //ファイルをクローズ．
fclose(outputdata_xyz);
return 0;
}
//バグ．複素数の負値の扱いがおかしい．
//また，mを負値にするとエラーで止まる．
//だけど，それっぽいグラフは書けた．
上の値設定の状態で，出力されるグラフは次のようになります．
早速ですが，上の関数が不完全であることがわかります．l=3, m=0のとき，
この関数はグラフの形がもう一つ上下逆に組み合わさった形になるはずなのですが，
それが実現できていないのです．何で？
見ず知らずの人に尋ねるのも何とも言えないのですが，わかる方，よろしければコメントくださいませんでしょうか．お待ちしております．
ギャラリーとして，他の条件のグラフも載っけときます．
それっぽくはなってますね……．
