oslo渐晕分析源代码(二)

#define PRECISION 1e-9

static int
vigchk(double max_pupil)
// hlp:  <P class="Edition">OSLO Premium/OSLO Standard</P>
// hlp:  <P>Routine to find <a href="../evaluate/illumination/vignetting.htm">fractional vignetting factors</a>.</P>
// hlp:  <P>Information is displayed for the <a href="../evaluate/setup/objectpoint.htm">current object point</a>.</P>
// hlp:  <P>The computed fxmax value is approximate, assuming that the greatest FX occurs
// hlp:  at the average FY point, i.e. (fymax + fymin)/2.</P>
// hlp:  <P>The algorithm code is available in CCL and can be modified for greater accuracy.</P>
// 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;
}