function [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans)

% ----[anton.public]
%
% function 
%
% [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans);
% 
%	inputs: 
%
%  i_detec:	=2: Xtomo prototype cameras (shot# < 6768)
%		=1: Xtomo 9-cameras	    (shot# > 682x)
%
%	outputs:
%
%  fans:    camera switch, 1=on,0=off (1x10)
%  vangle:  angle between detect. surface normal and pos. x-axis (1x10)
%  xchord:  two x-coordinates (2xnl) and
%  ychord:  two y-coord. for each line (2xnl), they specify start + end points
%  aomega:  etendue in mm^2 x steradians
%  angfact: angular factors, inverse of relative etendue (throughput) (20x10)
%
%	uses:
%		AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);
%		angular_fact_*.mat , '*'=i_detec
%  
%
%---------------- M.Anton 14/3/95 -------------------------------------------

disp('*----------------------------*')
disp('|   this is xtomo_geometry   |')
disp('*----------------------------*')

global xap yap xdet ydet
global ae da
% ======== tokamak parameters ================================================

load tcv_vesc1




xcont=rzvin(:,1);
ycont=rzvin(:,2);
xmin=min(xcont);
ymin=min(ycont);
xmax=max(xcont);
ymax=max(ycont);
xedge=100;
yedge=60;
          



% =========  detector parameters =============================================

if i_detec==2

	cw=1;                         % detector numbers cw=1:clockwise cw=0:ccw

	if nargin<2
		fans=[0 0 0 0 0 0 1 0 1 0];   % camera switch
	end

	vangle=[90 90 90 0 0 0 0 -90 -90 -90];
                                % angle of detector surface normal 
	xpos=[0 0 0 0 0 0  118.05 0 87.84 0]; 
                                % x position of the diaphragmas in [cm]
	ypos=[0 0 0 0 0 0 -46 0 -80.45 0];
                                % y position of the diaphragmas in [cm]
	ae=[0 0 0 0 0 0 -2.5 0 -0.1 0]/10; 
				%excentricity array/diaphragm in [cm]
	da=[0 0 0 0 0 0 10.1  0 -11.7 0]/10;   
                                 % diaphragma-array distance in [cm]
	da2=da;

	d1=0.950;                     % detector width      in mm
	d2=4.050;                     % detector length     in mm
	b1=1.000;                     % aperture width      in mm
	b2=4.000*ones(1,10);          % aperture length     in mm
	b3x=0;                        % aperture thickness  in mm
	b3y=0;

elseif i_detec==1

   cw=1;                         % detector numbers cw=1:clockwise cw=0:ccw

   if nargin<2	
   	fans=[1 1 1 1 1 1 1 1 0 1];   % camera switch
   end

   vangle=[90 90 90 0 0 0 0 -90 -90 -90];
                                 % angle of detector surface normal 
   xpos=[71.5 87.7 103.9 123.1  123.1 123.1 123.1 104.04 87.84 71.64]; 
                                 % x position of the diaphragmas in [cm]
   ypos=[80.35 80.35 80.35 48.1 1.95 -2.45 -48.6 -80.35 -80.35 -80.35];
                                 % y position of the diaphragmas in [cm]
   ae=[-8 0 8 5 9 -9 -5 8 0 -8]/10; %excentricity array/diaphragm [cm]


  ae= ae + [-0.0915 0   0.1361    0.2123    0.0923   -0.0994   ...
 		 0.0872   -0.1520  0  0.9410 ]/10;
  ae(1)=ae(1)+0.1/10;
  ae(3)=8/10+0.14/10;
  ae(4)=4.9/10;
  ae(5)=9/10+0.2/10;
  ae(6)=ae(6)-0.2/10;
  ae(7)=-4.9/10;
  ae(10)=-7.1/10;

   da=   [12.4 9.9 12.4 9.9 13.4 13.4  9.9   -12.4  -9.9 -12.4]/10;
                                 % diaphragma-array distance in [cm] (poloidal)
   da2=[37 34.4 37 55.9 59.4 59.4 55.9 37 34.4 37]/10;
                                 % dist to diaphragm in toroidal direction [cm];
   deltada=[   -0.0311 0  -0.0458   -0.1179   -0.0615   -0.1105 ...
			  -0.0510   -0.0515  0  -0.3223]/10;
  deltada(4)=0;
  deltada(6)=0;
                                                                 

   da=da+ deltada;

   da2=da2+deltada;


   d1=0.90;			 % detector width      in mm
   d2=4.0;			 % detector length     in mm
   b1=0.800;			 % aperture width      in mm (pol.)
   b2=[8 8 8 15 15 15 15 8 8 8];			 
				 % aperture length     in mm (tor.)
   b3x=0.020;			 % aperture thickness  in mm (poloidal)
   b3y=0;			 % aperture thickness  in mm (toroidal)
end



%======== calculation of the chords of view ===================================

nact=sum(fans);
iact=find(fans);
ndet=20;
ncam=10;



% ---- apertures: ------------------

xap=ones(ndet,1)*xpos(iact);
xap=xap(:)';
yap=ones(ndet,1)*ypos(iact);
yap=yap(:)';

% ---- detectors: ------------------

vorz(find(vangle==90))=(-1)^(cw+1)*ones(size(find(vangle==90)));
vorz(find(vangle==0))=(-1)^cw*ones(size(find(vangle==0)));
vorz(find(vangle==-90))=(-1)^cw*ones(size(find(vangle==-90)));


dete=(-9.025:0.950:9.025)'/10*vorz(iact)+ones(ndet,1)*ae(iact);

dum_ae=dete(:)';

dum_vangle=ones(ndet,1)*vangle(iact);
dum_vangle=dum_vangle(:)';



ivert=find(dum_vangle==90 | dum_vangle==-90);
ihori=find(dum_vangle==0);

dum_da=ones(ndet,1)*da(iact);
dum_da=dum_da(:)';

dxd=zeros(1,ndet*nact);
dyd=zeros(1,ndet*nact);

dxd(ivert)=dum_ae(ivert);
dxd(ihori)=dum_da(ihori);

dyd(ivert)=dum_da(ivert);
dyd(ihori)=dum_ae(ihori);

xdet=xap+dxd;
ydet=yap+dyd;


%plot_vessel(rzvin,rzvout)
%hold on
% plot(xap,yap,'.g',xdet,ydet,'.m')


% ---- calculate the equations of lines of sight

m=(ydet-yap)./(xdet-xap);
b=ydet-m.*xdet;

nl=length(xdet);
xchord=zeros(2,nl);
ychord=zeros(2,nl);


xchord(1,:)=xdet;ychord(1,:)=ydet;

iup=find(dum_vangle==90);
isi=find(dum_vangle==0);
ido=find(dum_vangle==-90);

if ~isempty(iup)
ychord(2,iup)=ymin*ones(size(iup));
xchord(2,iup)=(ychord(2,iup)-b(iup))./m(iup);
end
if ~isempty(ido)
ychord(2,ido)=ymax*ones(size(ido));
xchord(2,ido)=(ychord(2,ido)-b(ido))./m(ido);
end
if ~isempty(isi)
xchord(2,isi)=xmin*ones(size(isi));
ychord(2,isi)=m(isi).*xchord(2,isi)+b(isi);
end

ileft=find(xchord(2,:)<xmin);

if ~isempty(ileft)
xchord(2,ileft)=xmin*ones(size(ileft));
ychord(2,ileft)=m(ileft).*xchord(2,ileft)+b(ileft);
end
irig=find(xchord(2,:)>xmax);
if ~isempty(irig)
xchord(2,irig)=xmax*ones(size(irig));
ychord(2,irig)=m(irig).*xchord(2,irig)+b(irig);
end



%======== prepare output ======================================================

vangle=vangle(iact);

%======== calculation of angular correction factors, if necessary =============

if i_detec==2 & exist('angular_fact_2.mat')==2

		disp('loading angular_fact_2')
		load angular_fact_2


elseif i_detec==1 & exist('angular_fact_1.mat')==2

		disp('loading angular_fact_1')
		load angular_fact_1

else

		aomega=zeros(ndet,ncam);
		angfact=ones(ndet,ncam);


		for l=1:sum(fans)

%		Z0X=abs(da(iact(l))*10)
%	 	Z0Y=abs(da2(iact(l))*10)
%		X0=ae(iact(l))*10 % back to mm, sorry about that...
%		X0=X0*vorz(iact(l))
%		B2=b2(iact(l))
%		AOMEGA=etendue_n(b1,B2,b3x,b3y,Z0X,Z0Y,X0,cw);

		b1x=0.8;
		b1y=6;
		b1z=0.02;
		b2x=10000000;
		b2y=b2(iact(l));
		b2z=0;
		z01=abs(da(iact(l))*10);
		z02=abs(da2(iact(l))*10);
		X0=ae(iact(l))*10*vorz(iact(l));
keyboard
		AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);

			aomega(:,iact(l))=AOMEGA(:,1);

		
		end


		indm=min(find(aomega==max(aomega(:))));
		aomegan=aomega/aomega(indm);
		nonz=find(aomega);
		angfact(nonz)=ones(size(nonz))./aomegan(nonz);

		angfact=round(1000*angfact)/1000;


		unitstring='units aomega: mm^2 * sterad';

		if i_detec==1
			save angular_fact_1 angfact aomega unitstring
		elseif i_detec==2
			save angular_fact_2 angfact aomega unitstring
		end

end
return

th=atan(diff(ychord)./diff(xchord));
thp=th;
neg=find(thp<0);
thp(neg)=180+thp(neg);
thdet=ones(1,20*sum(fans));
for k=1:sum(fans)
     thdet((k-1)*20+1:k*20)=vangle(k)*ones(1,20);
end
angles=thdet-thp;
mist=find(angles<0 & abs(angles)>90);
angles(mist)=angles(mist)+180;
th_inc=angles*pi/180;


% ---- correct for the edges of tcv ( some chords may be too long )





	down=find(xcont>xedge & ycont<-yedge);
	up=find(xcont>xedge & ycont>yedge);
	cd=polyfit(xcont(down),ycont(down),1);
	cu=polyfit(xcont(up),ycont(up),1);


	iu1=find(xchord(1,:)>xedge & ychord(1,:)>0 & dum_vangle==-90 );
	if ~isempty(iu1)
		xchord(1,iu1)=-(b(iu1)-cu(2))./(m(iu1)-cu(1)+eps);
		ychord(1,iu1)=m(iu1).*xchord(1,iu1)+b(iu1);
	end

	iu2=find(xchord(2,:)>xedge & ychord(2,:)>0 & ychord(1,:) & ....
			dum_vangle==-90);
	if ~isempty(iu2)
		xchord(2,iu2)=-(b(iu2)-cu(2))./(m(iu2)-cu(1)+eps);
		ychord(2,iu2)=m(iu2).*xchord(2,iu2)+b(iu2);
	end

	id1=find(xchord(1,:)>xedge & ychord(1,:)<0 & dum_vangle==90);
	if ~isempty(id1)
		xchord(1,id1)=-(b(id1)-cd(2))./(m(id1)-cd(1)+eps);
		ychord(1,id1)=m(id1).*xchord(1,id1)+b(id1);
	end

	id2=find(xchord(2,:)>xedge & ychord(2,:)<0 & dum_vangle==90);
	if ~isempty(id2)
		xchord(2,id2)=-(b(id2)-cd(2))./(m(id2)-cd(1)+eps);
		ychord(2,id2)=m(id2).*xchord(2,id2)+b(id2);
	end

	ilow=find(ychord(1,:)<ymin);
	ihig=find(ychord(1,:)>ymax);
	ilef=find(xchord(1,:)<xmin);
	irig=find(xchord(1,:)>xmax);
	if ~isempty(ilow)
		ychord(1,ilow)=ymin*ones(size(ilow));
		xchord(1,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
	end
	if ~isempty(ihig)
		ychord(1,ihig)=ymax*ones(size(ihig));
		xchord(1,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
	end
	if ~isempty(ilef)
		xchord(1,ilef)=xmin*ones(size(ilef));
		ychord(1,ilef)=m(ilef)*xmin+b(ilef);
	end
	if ~isempty(irig)
		xchord(1,irig)=xmax*ones(size(irig));
		ychord(1,irig)=m(irig)*xmax+b(irig);
	end


	ilow=find(ychord(2,:)<ymin);
	ihig=find(ychord(2,:)>ymax);
	ilef=find(xchord(2,:)<xmin);
	irig=find(xchord(2,:)>xmax);
	if ~isempty(ilow)
		ychord(2,ilow)=ymin*ones(size(ilow));
		xchord(2,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
	end
	if ~isempty(ihig)
		ychord(2,ihig)=ymax*ones(size(ihig));
		xchord(2,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
	end
	if ~isempty(ilef)
		xchord(2,ilef)=xmin*ones(size(ilef));
		ychord(2,ilef)=m(ilef)*xmin+b(ilef);
	end
	if ~isempty(irig)
		xchord(2,irig)=xmax*ones(size(irig));
		ychord(2,irig)=m(irig)*xmax+b(irig);
	end





