const
     FaradayVoltage=550;
	 FaradayCount=0.642391669;	// SEC=0.642391669 @ 550V, SEC=120.721528 @ 850V
	 Intercept=-3.967736993;
	 Slope=0.006956141;

var
    Exposure,Voltage,Gain,Offset;
	Left,Right,Top,Bottom,Region;
	ROI,Size_x,Size_y;
	Img,RegImg,ImgName;
	SX,SY,TY;
	wX,wY;
	cX,cY,p;
	ADU,SI,Ne;
	pX,pY;

Function FWHM(A,npix);
		  var
		  	 i, iMax, Max, Min, Half, x1, x2, w;
		begin
			 Max:=A[0];
			 Min:=A[0];
			 x1:=0;
			 x2:=npix-1;
			 for i:=0 to npix-1 do
			 begin
			 	  if (Max<A[i]) then
				  begin 
				  	 Max:=A[i];
					 iMax:=i;
				  end;
				  if (Min>A[i]) then Min:=A[i];
			 end;
			 Half:=0.5*(Max+Min);
			 for i:=0 to npix-2 do
			 begin
			 	  if ((A[i]<=Half) and (Half<=A[i+1])) then x1:=i+(Half-A[i])/(A[i+1]-A[i]);
			 	  if ((A[i]>=Half) and (Half>=A[i+1])) then x2:=i+(Half-A[i])/(A[i+1]-A[i]);
			 end;
			 //WriteInfo('FWHM = ',x2-x1,' = ',x2,' - ', x1);
			 FWHM:=x2-x1;
		end;

Function AreaCenter(A,npix);
		  var
		  	 i, Min, sumy,sumxy, xc ;
		begin
			 sumy:=0;
			 sumxy:=0;
			 Min:=A[0];
			 for i:=0 to npix-1 do
			 begin
				  if (Min>A[i]) then Min:=A[i];
			 end;
			 for i:=0 to npix-1 do
			 begin
			 	  sumy:=sumy+(A[i]-Min);
				  sumxy:=sumxy+i*(A[i]-Min);
			 end;
			 //WriteInfo('FWHM = ',x2-x1,' = ',x2,' - ', x1);
			 AreaCenter:=sumxy/sumy;
		end;		

begin
     Offset:=28;
     Voltage:=550;
	 Exposure:=10;
//	 Exposure:=1;
	 SelectImage( 'Select IMAGE to process',Img ) ;
	 if IsNull( Img ) then halt('No image was selected') ;
	 ImgName:=GetName(Img);
	 ROI:=GetROI(Img);   // image and ROI starts from 0,0 on top left
	 if IsNull( ROI ) then 
	 begin
		GetXYSize(Img,Right,Bottom);
	    Left:=0;
		Top:=0;
		Right:=Right-1;
		Bottom:=Bottom-1;
		WriteStatus('No ROI is defined in '+ImgName);
     end
	 else
	 begin
        Left:=ROI[0,0];
        Top:=ROI[1,0];
        Right:=ROI[2,0];
        Bottom:=ROI[3,0];
     end;
	 if (Right-Left)<3 then WriteInfo('The left two columns will be cut ...');
		
	 WriteStatus('ROI is '+Str(Left)+','+Str(Right)+', '+Str(Top)+','+Str(Bottom));
	 GetNumber('Enter CCD offset ADU',Offset);
	 GetNumber('Enter ImageIntensifier setting voltage',Voltage);
	 GetNumber('Enter Exposure time in millisecond',Exposure);
	 Gain:=FaradayCount*Exp(Slope*(Voltage-FaradayVoltage)*LN(10));
	 //RegImg:=GetRegion(Img,Left-1,Top-1,Right-1,Bottom-1);		  // obsolete function
	 RegImg:=CreateImage(typ_Double,Right-Left+1-2,Bottom-Top+1);
	 RegImg[0..(Right-Left-2),0..(Bottom-Top)]:=Img[(Left+2)..(Right),(Top)..(Bottom)]; // +2 eliminates CCD artifact
	 Size_x:=GetXSize(RegImg);
	 Size_y:=GetYSize(RegImg);
	 SY:=MeanOverX(RegImg);
	 SX:=MeanOverY(RegImg);
	 SI:=SumOf(RegImg);
	 Ne:=(SI-Offset*Size_x*Size_y)/(Gain*Exposure);
	 ADU:=MeanOf(RegImg);

	 //WriteInfo('Size of ' + ImgName + ' = ' , Size_x , ' x ' , Size_y);

 	 wX:=FWHM(SX,Size_x);
	 cX:=2+Left+AreaCenter(SX,Size_x);
	 pX:=CreatePlot('Row Profile');
	 Plot(pX,SX);
	 
     TY:=CreateArray(typ_Double,Size_y,1);
	 TY[..]:=SY[0,..];
	 wY:=FWHM(TY,Size_y);
	 cY:=Top+AreaCenter(TY,Size_y);
	 pY:=CreatePlot('Col Profile');
	 Plot(pY,TY);
	 
     WriteStatus('#e = ',Ne,' , FWHM = ',wX,' x ',wY,' @ ',cX,' , ',cY);
	 WriteInfo('#e = ',Ne,' , FWHM = ',wX,' x ',wY,' @ ',cX,' , ',cY);
	 
	 Free(Img);
	 Free(RegImg);
	 Free(ROI);
	 Free(TY);
	 Free(SY);
	 Free(SX);
end
