oslo渐晕分析源代码(二)

╰+哭是因爲堅強的太久メ 2022-05-11 07:20 178阅读 0赞

#define PRECISION 1e-9

static int
 vigchk(double max_pupil)
 // hlp: OSLO Premium/OSLO Standard
 // hlp: Routine to find <a href="../evaluate/illumination/vignetting.htm">fractional vignetting factors</a>.
 // hlp: Information is displayed for the <a href="../evaluate/setup/objectpoint.htm">current object point</a>.
 // hlp: The computed fxmax value is approximate, assuming that the greatest FX occurs
 // hlp: at the average FY point, i.e. (fymax + fymin)/2.
 // hlp: The algorithm code is available in CCL and can be modified for greater accuracy.
 // kwd: vignetting, apertures, field points
 // cat: analysis tools
 {
 	double fymin, fymax, fx, fytmp;
 	double _fby, _fbx;
 	double fymin_sup, fymin_inf, fymax_sup, fymax_inf, fx_sup, fx_inf;
 	double fymin_sup_sav, fymax_inf_sav;
 	char found_ray_through;
 	int nbr_rays;
 	int ssb_row_sav;
 	int i, j, k;
 	
 	SAVE_DISPLAY_PREFS;
 	fymin = -max_pupil;
 	fymax = +max_pupil;
 	
 	stp outp off;
 	ssb_row_sav = sbrow();
 	
 	fytmp = 0;
 	
 	// Find a meridional ray that gets through system
 	// 下面是将孔径坐标-1到1，先分割成20份，从(-1，0)追迹到(1,0)
 	// 直到发现第一条可穿过系统的子午光线，fytmp存下这个第一条光线的孔径坐标
 	nbr_rays = 10;
 	found_ray_through = FALSE;
 	while (!found_ray_through && 1/nbr_rays > PRECISION/10)
 	{
 		for (i = -NBR_RAYS; i < NBR_RAYS; i++)
 		{
 			if (!found_ray_through)
 			{
 				//将即将输入的spreadbuffer的两行，填充为1e20
 				ssbuf_reset(ssb_row_sav, 2, 1e20);
 				trace_ray i/NBR_RAYS 0.0;
 				//如果这条光线能穿过系统的话，ssb(1,1)存的是光线在像面上交点的y坐标
 				//若无法穿过系统的话，将什么也不输出
 				if (ssb(1, 1) != 1e20)
 				{
 					found_ray_through = TRUE;
 					fytmp = i/nbr_rays;
 				}
 				ssbuf_reset(-ssb_row_sav, 0);
 			}
 			else
 				break;
 		}
 		//如果没发现这条穿过系统的子午光线，则将光阑多划分十倍
 		//重复以上过程，继续寻找，直到找到
 		//或划分的格子已经小于1e-9，实在是找不到了，停止搜寻
 		nbr_rays *= 10;
 	}
 	
 	if (!found_ray_through)
 	{
 		aprintf("Could not get a ray through system\n");
 		RESTORE_DISPLAY_PREFS;
 		return NO_RAY_THROUGH;
 	}
 	
 	// Find fymin
 	// fymin_sup存刚刚找到的第一条可穿过系统的子午光线
 	// 用二分法，在[-max_pupil,fytmp]寻找能穿过系统的光线的最小孔径坐标，存到fymin中
 	fymin_sup = fytmp;
 	fymin_inf = -max_pupil;
 	while (fymin_sup - fymin_inf > PRECISION)
 	{
 		ssbuf_reset(ssb_row_sav, 2, 1e20);
 		trace_ray((fymin_inf+fymin_sup)/2, 0.0);
 		if (ssb(1, 1) == 1e20)
 			fymin_inf = (fymin_inf + fymin_sup)/2;
 		else
 		{
 			fymin_sup_sav = fymin_sup;
 			fymin_sup = (fymin_inf + fymin_sup)/2;
 		}
 		ssbuf_reset(-ssb_row_sav, 0);
 	}
 	fymin = fymin_sup_sav;
 	
 	// Find fymax
 	// 用二分法，在[fytmep,max_pupil]寻找能穿过系统的子午光线的最大孔径坐标，存到fymax中
 	fymax_sup = +max_pupil;
 	fymax_inf = fytmp;
 	while (fymax_sup - fymax_inf > PRECISION)
 	{
 		ssbuf_reset(ssb_row_sav, 2, 1e20);
 		trace_ray((fymax_inf+fymax_sup)/2, 0.0);
 		if (ssb(1, 1) == 1e20)
 			fymax_sup = (fymax_inf + fymax_sup)/2;
 		else
 		{
 			fymax_inf_sav = fymax_inf;
 			fymax_inf = (fymax_inf + fymax_sup)/2;
 		}
 		ssbuf_reset(-ssb_row_sav, 0);
 	}
 	fymax = fymax_inf_sav;
 	
 	
 	// Find the fractional aperture fxmin of the first sagital ray @ (fymin+fymax)/2
 	// 在子午坐标为(fymin+fymax)/2处，搜寻能穿过系统的弧矢光线的最大弧矢坐标
 	// 用二分法，在[0,max_pupil]之间找
 	fx_sup = +max_pupil;
 	fx_inf = 0;
 	while (fx_sup - fx_inf > PRECISION)
 	{
 		ssbuf_reset(ssb_row_sav, 2, 1e20);
 		trace_ray((fymin+fymax)/2, (fx_inf+fx_sup)/2);
 		if (ssb(1, 1) == 1e20)
 			fx_sup = (fx_inf + fx_sup)/2;
 		else
 			fx_inf = (fx_inf + fx_sup)/2;
 		ssbuf_reset(-ssb_row_sav, 0);
 	}
 	fx = (fx_inf + fx_sup)/2;
 	
 	ssbuf_reset(ssb_row_sav, 3);
 	stp outp off;
 	trr;
 	_fby = ssb(1, 1);
 	_fbx = ssb(1, 2);
 	ssbuf_reset(-ssb_row_sav, 3);
 	RESTORE_DISPLAY_PREFS;
 	aprintf("\n@(vigchk)*VIGCHK CFG FBY FBX FYMIN FYMAX APPROX_FXMAX\n");
 	aprintf(" %d %.2f %.2f % .3f % .3f % .3f \n", cfg, _fby, _fbx, fymin, fymax, fx);
 	//如果含有特殊孔径的，可能不准确
 	k = 0;
 	for (j = 1; j <= ims; j++)
 		k += numsap[j];
 	if (k)
 	{
 		aprintf("Warning: This system contains special aperture data.\n");
 		aprintf(" Vignetting may not be correct.\n");
 	}
 	
 	return 0;
 }