光学设计的秘密花园

　　张大千临摹了几百张明代石涛的山水，画出的画以假乱真，后来他去敦煌潜心临摹几年，回来画风大变，终成大家。

　　那我就转载一下百度百科的词条警醒一下自己。

...

　　对于石涛他尤为推崇，他用了大量精力去学习石涛的绘画艺术，把石涛的艺术比喻成万里长城。他还从李瑞清之弟李筠庵学会仿制古画的方法，并做了许多石涛的赝品，曾多次骗过程霖生。石涛画境变化无尽，新颖怪奇而又法度严谨，大千正是通过石涛而涵泳了唐宋元明百家之长。在二十年代的上海他获得了"石涛专家"的美誉。

...

　　40年代，张大千自费赴敦煌，耗时三年大量临摹了石窟壁画，并将之宣传介绍，使敦煌艺术宝库从此为国人和世界广为瞩目。从此，张大千的画风也为之一变，善用复笔重色，高雅华丽，潇洒磅礴，被誉为“画中李白”、“今日中国之画仙”。1942年，春末，他决定举家赴敦煌临摹壁画。在敦煌的生活是艰苦的，恶劣的气候条件，再加上住在与世隔绝的石洞子里对于一般人是无法忍受的，但对于一个艺术家却可能是有益的。大千在此时画风为之一变，他善用复笔重色，笔力也变得丰厚浓重。可以说在敦煌时期对他风格形成起到了至关重要的作用。他不仅考察莫高窟的壁画佛像，还对莫高窟进行了编号，成为了为莫高窟编号的第一人，为保存文化遗产做了积极的贡献。他还到达青海西宁，邀请藏族画师共同赴敦煌协助临摹工作。为尊重他人劳动成果和对摹品的负责，每幅画上都注明了画家的名字，凡与他合作也注明了作者的名字，所以在他许多临摹的敦煌壁画中都标有"番僧某某同画"。两年后他结束了在莫高窟的临摹工作，同时开始了对安西榆林窟的临摹工作。

　　在离开莫高窟后，他花了大量时间对十六国、北魏、北周、隋、唐、五代、宋、西夏、元各朝的壁画代表作及雕塑进行了临摹，共有摹品共记三百多幅。...

　　从入行以来，一直都有这样几个疑问，盘桓在我的心头，挥之不去。

　　一， 高斯光束通过实际透镜的情况下，高斯光束变换和像差的综合作用是什么？

　　二， 厄米高斯光束和拉盖儿高斯光束在用几何光学的处理方法上有什么不同？

　　三， 基于光线追迹的变换处理方法得到物理光学演算结论的能力极限在哪里？

　　在这里，我刻下一个对明年的承诺：明年年底的时候，我要解决这三个疑问。

　　赤手缚虎？不，我不是周处。我需要找到一些工具。首先就是“the ray and wave theory of lenses”。最早知道这本书，是陈志隆老师的推荐。

…报告完这本书之后, 会再念另一本书；念的是A. Walther, "the ray and wave theory of lenses",(习题要作;写程序的要用Mathematica与OSLO's ccl或dll programming);之后才会读期刊论文.谈作论文研究.

--- Walther的书是作进阶开发的基础

--- Mathematica 是为了symbolic calculation用

--- OSLO's ccl 与 dll programming 是为作 optical module 研发的implementation用

…

　　看起来这本书很是被重视。究竟内容如何呢，我在这里下载了部分内容，并翻译在下面。

　　现代光学的很多书籍，常常局限于近轴范围来处理透镜理论。像差则被当作近轴理想的一个补丁。这些书中一直没有解释清楚的是，像差的存在是由于物理法则还是由于设计者有限的技艺。读者一直被灌输这样的观点，镜头的设计必须尽可能的遵循近轴光学的成像规律。这种说法令人遗憾。因为从十八世纪我们就知道，如果将近轴光学的公式用于有限的高度和角度，则会导出投影几何。这是数学中一个重要分支。但对于实际透镜的行为而言，这是一个糟糕的描述。

　　过去四十年傅立叶光学的发展让透镜理论和物理实际之间更加接近。但由于许多作者不顾及实验室中大孔径镜头的普遍存在，在处理成像公式理论时满足于近轴近似，导致有很多光学见解的遗漏。笔者认为，要么能够提供一个清晰且令人信服的验证，说明近轴近似可靠的适用于大孔径和大视场；要么应该可靠的发展一套没有近轴近似的光学理论。这套可靠的理论是存在的，不过被埋没于充满细节描述的书籍和论文中，以致于初学者迷失于数学描述的错综复杂之中。

　　本书的主要目的是描述古老而完备的eikonal函数可以用来发展一套镜头的理论，可同时在几何光学和波动光学成像理论中，很好的处理于大孔径大视场的情况。笔者尽力遵循一下三条原则：

• 几何光学和物理光学之间的联系不能只是无实际意义的数学推导，而必须是一个有力的工程工具，使得读者可以在几何光学和物理光学中互相推导切换。
• 能量守恒必须是这套理论的组成部分。避免出现神秘而不确定的振幅因子。
• 不考虑只是优雅的理论。任何获得的结果必须有助于讨论更深入，或者有其实际应用。

　　很好，就是它了！

　　《透镜设计基础》中第314页到第321页详述了用手工计算的方法设计一个Cooke式三片镜，也就是K. Schwarzschild方法。在以下的内容中，首先演示如何cooke_triplet命令进行Cooke式三片镜的设计，然后讲解如何将这一方法整合到OSLO的环境中。

　　启动OSLO，在lens spreadsheet关闭时按Ctrl+D，出现如下对话框，每一参数的意义请参照《透镜设计基础》中的解释。

　　点选OK，然后在命令行中依照提示依次输入三个镜片的中心厚度，比如0.4，0.25，0.45，分别按回车键确认。这时，当前的OSLO环境中已经生成了一个符合参数要求的Cooke式三片镜，分析图如下。

　　从图中可以看出子午像面过于前倾，可以适当调整一下。将AC的预设值从-0.09调整到-0.06。

　　分析新的镜头，可以看到子午像面不再严重前倾，像散有所改善。

　　讲解如何将这一方法整合到OSLO的环境中。首先打开\private\ccl\a_menu.ccl文件，加入下图中的高亮部分。

　　再打开\private\ccl\a_global.ccl加入下列内容。

double y1 {default(noquery) = 1.111111};

double y2 {default(noquery) = 0.888888};

double y3 {default(noquery) = 0.999999};

double rn1 {default(noquery) = 1.620318};

double vn1 {default(noquery) = 60.31};

double rn2 {default(noquery) = 1.616445};

double vn2 {default(noquery) = 36.62};

double rn3 {default(noquery) = 1.620318};

double vn3 {default(noquery) = 60.31};

double total_efl {default(noquery) = 10};

double lch {default(noquery) = -0.02};

double Petzval {default(noquery) = 0.035};

double K_aim {default(noquery) = 1};

double fie_ang {default(noquery) = 20};

int step {default(noquery) = 500};

double c1_start {default(noquery) = 0.2};

double c1_end {default(noquery) = 0.4};

double SC_target {default(noquery) = -0.08};

double CC_target {default(noquery) = 0.0025};

double AC_target {default(noquery) = -0.09};

double TchC_target {default(noquery) = 0.0};

double Wav1 {default(noquery) = 0.5892938};

double Wav2 {default(noquery) = 0.4861327};

double Wav3 {default(noquery) = 0.6562725};

最后是函数的主体部分，编译后就可以了。

cmd scalecv(int start_surf, int end_surf, float efl)

{

float feedback=1.0;

stp outp off;

int counter;

do{

for(counter=start_surf;counter<end_surf+1;counter++)

cv(counter,cv[counter]*feedback);

feedback=pe_efl(start_surf,end_surf)/efl;

}

while(fabs(feedback-1)>0.00000001);

stp outp on;

}

cmd Gsum ( float refra_index, float result[] )

{

result[0] = 0;

result[1] = refra_index * refra_index * ( refra_index - 1 ) / 2 ;

result[2] = ( 2 * refra_index + 1 ) * ( refra_index - 1 ) / 2 ;

result[3] = ( 3 * refra_index + 1 ) * ( refra_index - 1 ) / 2 ;

result[4] = ( refra_index + 2 ) * ( refra_index - 1 ) / ( 2 * refra_index ) ;

result[5] = 2 * ( refra_index + 1 ) * ( refra_index - 1 ) / refra_index ;

result[6] = ( 3 * refra_index + 2 ) * ( refra_index - 1 ) / ( 2 * refra_index ) ;

result[7] = ( 2 * refra_index + 1 ) * ( refra_index - 1 ) / ( 2 * refra_index );

result[8] = refra_index * ( refra_index - 1 ) / 2;

}

double SC_coef[3];

double CC_coef[3];

double AC_coef[3];

double SC_ap_coef[3];

double CC_ap_coef[3];

double AC_ap_coef[3];

cmd sc_sei(double sc_sei_coef[], double y_para, double u_prime_0, double G_sum[], double curvature, double inv_objdist)

{

//calculate the constant

double temp_coef_for_all;

int ii;

temp_coef_for_all = -1 * y_para * y_para * y_para * y_para / ( u_prime_0 * u_prime_0 );

sc_sei_coef[0] = G_sum[1] * curvature * curvature * curvature;

sc_sei_coef[1] = -1 * G_sum[2] * curvature * curvature;

sc_sei_coef[0] += G_sum[3] * curvature * curvature * inv_objdist;

sc_sei_coef[2] = G_sum[4] * curvature;

sc_sei_coef[1] += -1 * G_sum[5] * curvature * inv_objdist;

sc_sei_coef[0] += G_sum[6] * curvature * inv_objdist * inv_objdist;

for (ii = 0; ii < 3; ii++)

{

sc_sei_coef[ii] *= temp_coef_for_all;

}

}

cmd cc_sei(double cc_sei_coef[], double y_para, double h_prime_0, double G_sum[], double curvature, double inv_objdist)

{

double temp_coef_for_all;

int ii;

temp_coef_for_all = -1 * y_para * y_para * h_prime_0;

cc_sei_coef[1] = 0.25 * G_sum[5] * curvature;

cc_sei_coef[0] = -1 * G_sum[7] * curvature * inv_objdist;

cc_sei_coef[0] += -1 * G_sum[8] * curvature * curvature;

cc_sei_coef[2] = 0;

for (ii = 0; ii < 3; ii++)

{

cc_sei_coef[ii] *= temp_coef_for_all;

}

}

cmd ac_sei(double ac_sei_coef[], double h_prime_0, double fl )

{

ac_sei_coef[0] = -0.5 * h_prime_0 * h_prime_0 / fl;

ac_sei_coef[1] = 0;

ac_sei_coef[2] = 0;

}

cmd sc_ap(double ap_shift[], double ap_attached[])

{

for (ii = 0; ii < 3; ii++)

{

ap_shift[ii] = ap_attached[ii];

}

}

cmd cc_ap(double cc_ap_shift[], double cc_ap_attached[], double sc_ap_attached[], double q_factor, double u_prime_0)

{

for (ii = 0; ii < 3; ii++)

{

cc_ap_shift[ii] = cc_ap_attached[ii] + sc_ap_attached[ii] * q_factor * u_prime_0;

}

}

cmd ac_ap(double ac_ap_shift[], double ac_ap_attached[], double cc_ap_attached[], double sc_ap_attached[], double q_factor, double u_prime_0)

{

for (ii = 0; ii < 3; ii++)

{

ac_ap_shift[ii] = ac_ap_attached[ii] + cc_ap_attached[ii] * 2 * q_factor / u_prime_0 + sc_ap_attached[ii] * q_factor * q_factor;

}

}

cmd cooke_triplet(double y1, double y2, double y3, double wav1, double wav2, double wav3, double rn1, double vn1, double rn2, double vn2,

double rn3, double vn3, double total_efl, double SC_target, double CC_target, double AC_target,

double Petzval, double lch, double TchC_target, double K_aim, double fie_ang, double c1_start, double c1_end, int step)

{

stp outp off;

double phi[3];

double y[3];

y[0] = y1; y[1] = y2; y[2] = y3;

double v[3];

v[0] = vn1; v[1] = vn2; v[2] = vn3;

double n[3];

n[0] = rn1; n[1] = rn2; n[2] = rn3;

double phi_total;

phi_total = 1/total_efl;

static double u0_obj = 0, u0_ims;

u0_ims = y1/total_efl;

double ptz;

ptz=Petzval;

double coef_phi[3][3], verctor_aber[3];

double ua, ub, uc;

ua=u0_obj;

double airspace_1, airspace_2;

double K;

double k_new, delta_yb, delta_yc;

double upr[4], ypr[3];

upr[0] = tan(-1*fie_ang*dr);

double TchC[3];

double TchC_total, TchC_total_new, TchC_aim;

TchC_aim = TchC_target;

for ( ii = 0, TchC_total = -0.02, delta_yc = 0.01, k = 0.5, delta_yb = 0.01 ; ii < 200; ii++)

{

y[2] += delta_yc;

for ( jj = 0; jj < 200; jj++)

{

y[1] += delta_yb;

for ( i = 0; i < 3; i++) coef_phi[0][i] = y[i];

for ( i = 0; i < 3; i++) coef_phi[1][i] = y[i] * y[i] / v[i];

for ( i = 0; i < 3; i++) coef_phi[2][i] = 1 / n[i];

verctor_aber[0] = y[0] * phi_total;

verctor_aber[1] = -1 * lch * u0_ims * u0_ims;

verctor_aber[2] = Ptz;

slveqs( coef_phi, phi, verctor_aber);

/////////////////////////

//ub uc is at object-side

/////////////////////////

ub = ua + y[0] * phi[0];

uc = ub + y[1] * phi[1];

airspace_1 = ( y[0] -y[1]) / ub;

airspace_2 = ( y[1] -y[2]) / uc;

K_new = airspace_1 / airspace_2;

if ( fabs( K_new - K_aim ) < 0.000001 ) break;

delta_yb = ( K_aim - K_new ) * delta_yb / ( K_new - K );

k = k_new;

}

ypr[0] = airspace_1 * upr[0] / ( 1 - airspace_1 * phi[0] ); ypr[1] = 0; ypr[2] = -1 * ypr[0] / k_new;

upr[1]= upr[0] + ypr[0] * phi[0];

upr[2]= upr[1] + ypr[1] * phi[1];

upr[3]= upr[2] + ypr[2] * phi[2];

TchC[0] = ( -1 * y[0] * ypr[0] * phi[0] ) / ( V[0] * u0_ims);

TchC[1] = 0;

TchC[2] = ( -1 * y[2] * ypr[2] * phi[2] ) / ( V[2] * u0_ims);

TchC_total_new = TchC[0] + TchC[1] + TchC[2];

if ( fabs( TchC_total_new - TchC_aim ) < 0.000001 ) break;

delta_yc = ( TchC_aim - TchC_total_new ) * delta_yc /( TchC_total_new - TchC_total );

TchC_total = TchC_total_new;

}

int lensID;

double g_coef[9];

double sc_coef[3][3];

double cc_coef[3][3];

double ac_coef[3][3];

double sc_ap_coef[3][3];

double cc_ap_coef[3][3];

double ac_ap_coef[3][3];

double temp_sc_coef[3];

double temp_cc_coef[3];

double temp_ac_coef[3];

double temp_sc_ap_coef[3];

double temp_cc_ap_coef[3];

double temp_ac_ap_coef[3];

double u0_=u0_ims;

double c[3];

c[0] = phi[0]/(n[0]-1);

c[1] = phi[1]/(n[1]-1);

c[2] = phi[2]/(n[2]-1);

double inv[3];

inv[0] = ua/y[0];

inv[1] = ub/y[1];

inv[2] = uc/y[2];

double h0_;

h0_= (ypr[2]-upr[3]*(y[2]/(uc+y[2]*phi[2])));

double q[3];

q[0] = ypr[0]/y[0];

q[1] = ypr[1]/y[1];

q[2] = ypr[2]/y[2];

for (lensID=0; lensID<3; lensID++)

{

gsum(n[lensID], g_coef);

//for(i=1; i<9; i++) prt g_coef[i];

//prt "=============================";

sc_sei(temp_sc_coef, y[lensID], u0_, g_coef, c[lensID], inv[lensID]);

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

{

sc_coef[i][lensID]=temp_sc_coef[i];

//aprt sc_coef[i][lensID];

};

//prt;

cc_sei(temp_cc_coef, y[lensID], h0_, g_coef, c[lensID], inv[lensID]);

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

{

cc_coef[i][lensID]=temp_cc_coef[i];

//aprt cc_coef[i][lensID];

};

//prt;

ac_sei(temp_ac_coef, h0_, (1/phi[lensID]));

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

{

ac_coef[i][lensID]=temp_ac_coef[i];

//aprt ac_coef[i][lensID];

};

//prt;

sc_ap(temp_sc_ap_coef, temp_sc_coef);

for(i=2; i>-1; i--)

{

sc_ap_coef[i][lensID]=temp_sc_ap_coef[i];

//aprintf("%f\t", sc_ap_coef[i][lensID]);

};

//prt;

cc_ap(temp_cc_ap_coef, temp_cc_coef, temp_sc_coef, q[lensID], u0_);

for(i=2; i>-1; i--)

{

cc_ap_coef[i][lensID]=temp_cc_ap_coef[i];

//aprintf("%f\t", cc_ap_coef[i][lensID]);

};

//prt;

ac_ap(temp_ac_ap_coef, temp_ac_coef, temp_cc_coef, temp_sc_coef, q[lensID], u0_);

for(i=2; i>-1; i--)

{

ac_ap_coef[i][lensID]=temp_ac_ap_coef[i];

//aprintf("%f\t", ac_ap_coef[i][lensID]);

};

//prt;

}

double surf_curv[6];

double AC_1_2;

double c5_coef[3];

double a,b,cc;

double c5[3];

double CC_1_3;

double c3;

double SC_total;

int srfID;

for (ii=1; ii<step; ii++)

{

surf_curv[0] = c1_start + ii*(c1_end - c1_start)/(step-1) ;

//prt surf_curv[0];

//change the C1, to get the C5;

//caculate AC*_1

AC_1_2 = ac_ap_coef[2][0] * surf_curv[0] * surf_curv[0] + ac_ap_coef[1][0] * surf_curv[0] + ac_ap_coef[0][0];

AC_1_2 += ac_ap_coef[0][1];

//prt "AC*_1_2=" AC_1_2;

c5_coef[2] = ac_ap_coef[2][2]; c5_coef[1] = ac_ap_coef[1][2];

c5_coef[0] = ac_ap_coef[0][2] + AC_1_2 - AC_target;

//prt "C5^2, C5^1, C5^0" (c5_coef[2]/c5_coef[2]) (c5_coef[1]/c5_coef[2]) (c5_coef[0]/c5_coef[2]);

a=c5_coef[2]; b=c5_coef[1]; cc=c5_coef[0];

//prt (b*b-4*a*cc);

if ( b*b-4*a*cc < 0 )

prt "No solution.";

else { c5[2] = (-1*b + sqrt(b*b-4*a*cc))/(2*a); c5[1] = (-1*b - sqrt(b*b-4*a*cc))/(2*a);}

//choose the less bending one, coef_calculating can be used

c5[0]=c5[1];

CC_1_3 = cc_ap_coef[2][0]*surf_curv[0]*surf_curv[0] + cc_ap_coef[1][0]*surf_curv[0] + cc_ap_coef[0][0];

CC_1_3 += cc_ap_coef[2][2]*c5[0]*c5[0] + cc_ap_coef[1][2]*c5[0] + cc_ap_coef[0][2];

//prt "CC_1_3=" CC_1_3;

c3 = (-1*cc_ap_coef[0][1]+(-1*CC_1_3)+CC_target)/cc_ap_coef[1][1];

//prt "C3=" c3;

sc_total = sc_ap_coef[2][0]*surf_curv[0]*surf_curv[0] + sc_ap_coef[1][0]*surf_curv[0] + sc_ap_coef[0][0];

sc_total += sc_ap_coef[2][1]*c3*c3 + sc_ap_coef[1][1]*c3 + sc_ap_coef[0][1];

sc_total += sc_ap_coef[2][2]*c5[0]*c5[0] + sc_ap_coef[1][2]*c5[0] + sc_ap_coef[0][2];

if (fabs(sc_total-SC_target)<0.001)

{

aprt "total SC = " sc_total;

aprt "C1, C3, C5" surf_curv[0] c3 c5[0];

break;

}

}

double th1;

double th2;

double th3;

input th1 "TH1";

input th2 "TH2";

input th3 "TH3";

len new "Cooke draft" total_efl;

ebr y1;

ang 1e-20;

des "ZhangJing";

uni 1.0;

sno1 "Produced by design_cooke.ccl";

wv wav1 wav2 wav3;

WW 1.0 1.0 1.0;

//srf 0

air;

ap 2e+19;

th 1.0e+20;

//srf 1

NXT;

gla mod "glass1" n[0] v[0];

cv surf_curv[0];

th th1;

ap noc y[0];

//srf 2

NXT;

air;

cv surf_curv[0]-c[0];

th airspace_1;

ap noc y[0];

//srf 3

NXT;

gla mod "glass3" n[1] v[1];

cv c3;

th th2;

ap noc y[1];

//srf 4

NXT;

air;

cv c3-c[1];

th airspace_2;

ap noc y[1];

ast;

//srf 5

NXT;

gla mod "glass5" n[2] v[2];

cv c5[0];

th th3;

ap noc y[2];

//srf 6

NXT;

air;

cv c5[0]-c[2];

th 0.0;

ap noc y[2];

//IMS

NXT;

air;

ap noc 0.0;

end;

stp outp on;

aprt (1/phi[0]);

aprt (1/phi[1]);

aprt (1/phi[2]);

scalecv(1,2, 1/phi[0]);

scalecv(3,4, 1/phi[1]);

scalecv(5,6, 1/phi[2]);

nodp(1,2,1);

e=ssb(sbrow-1,4);

nodp(3,4,1);

f=ssb(sbrow-1,3);

h=airspace_1+e-f;

th(2, h);

nodp(3,4,1);

e=ssb(sbrow-1,4);

nodp(5,6,1);

f=ssb(sbrow-1,3);

h=airspace_2+e-f;

th(4, h);

tra std loc All usr 1 0 n 1;

for (ii =1; ii<ims; ii++) ap(ii, chk, ssb(sbrow-8,1)+0.000001);

auf pxl;

ang fie_ang;

f 4 ins 1.2 0 0 0 0 -1.0 1.0 -1.0 1.0 1.0 0 0;

vig y y 5;

dlnr 1 0;

drw ims on;

stp outp on;

}