program glaremapmodel
	implicit none
	real*4 glaremodel(90,180),glaredge(90,180)
	real*4 tx,ty,cx,cy,cz,arg,c1,s1,c2,s2,s1wall,s2wall,cwall,cx1,cz1,a4,a5,a6,a7,a8,arg1,z4,z5,z6,z7,zz4,zz5,zz6	!,b
	integer*4 i,j,ii,jj,ibad
	data glaremodel /16200*0./
	data glaredge /16200*0./
	data s1wall,s2wall /60.,0./
c	data z4,z5,z6,z7,zz4,zz5,zz6  /1.2,.7,.5,.5,1.5,.9,.6/	!values 10 Oct 2005: too much Z6
c	data a8,a7,a6,a5,a4 /20., 0.04, 0.05, 0.1, 0./	!values 10 Oct 2005: too much Z6
	data z4,z5,z6,z7,zz4,zz5,zz6  /2.0,1.9,0.25,0.4,2.0,1.9,0.3/	!values 20 Oct 2005: Z6 adjusted down, Z4, Z5 adjusted up, b eliminated
	data a8,a7,a6,a5,a4 /20., 0.05, 0.025, 0.02, 0./	!13 Oct 2005	!note a8 is much bigger than the others: it alone offsets cz
c	data b /0.1/
	c1=cosd(s1wall)				!s1wall is the angle between baffle axis and normal plane to sept bottoms
	s1=sind(s1wall)
	c2=cosd(s2wall)				!s2wall is an ad hoc angle to move the hot spot in the + Theta-Y direction
	s2=sind(s2wall)
	ibad=0
c	
	do ii=1,90
	i=91-ii
	if(i.eq.86)then
	do j=10,171
		glaredge(i,j)=0.2			!Z7	!Mark verticals in edge map
	enddo
	endif
	if(i.eq.81)then
	do j=37,144
		glaredge(i,j)=0.6			!Z6
	enddo
	endif
	if(i.eq.70)then
	do j=60,121
		glaredge(i,j)=1.0			!Z5
	enddo
	endif
	if(i.eq.65)then
	do j=66,115
		glaredge(i,j)=1.4			!Z4
	enddo
	endif
	if(i.eq.58)then
	do j=72,109
		glaredge(i,j)=1.9			!Z3
	enddo
	endif
c	
	do j=91,180
c
	if(i.lt.86.and.j-90.eq.81)glaredge(i,j)=0.2	!Mark horizontals in edge map	
	if(i.lt.81.and.j-90.eq.54)glaredge(i,j)=0.6	
	if(i.lt.70.and.j-90.eq.31)glaredge(i,j)=1.0	
	if(i.lt.65.and.j-90.eq.25)glaredge(i,j)=1.4
	if(i.lt.58.and.j-90.eq.19)glaredge(i,j)=1.9
	if(i.lt.86.and.j-90.eq.81)glaredge(i,181-j)=0.2
	if(i.lt.81.and.j-90.eq.54)glaredge(i,181-j)=0.6	
	if(i.lt.70.and.j-90.eq.31)glaredge(i,181-j)=1.0	
	if(i.lt.65.and.j-90.eq.25)glaredge(i,181-j)=1.4
	if(i.lt.58.and.j-90.eq.19)glaredge(i,181-j)=1.9
c		
		jj=j-90
		tx=float(i)-0.5
		if(tx.eq.89.5)tx=89.
		ty=float(jj)-0.5
		if(ty.eq.89.5)ty=89.
		cx=sind(tx)/sqrt(1.+((tand(ty))**2) *(cosd(tx)**2))		!here find direction cosines
		cy=sind(ty)/sqrt(1.+((tand(tx))**2) *(cosd(ty)**2))
		arg=cx*cx+cy*cy
		cz=0.
		if(arg.ge.1.0)print *,i,j,arg
		if(arg.lt.1.0)cz=sqrt(1.-arg)
c		cwall=s1*cx+c1*cz						!find normal to septum-bottom planes
c		wallsave(i,j)=cwall
c		wallsave(i,181-j)=cwall						!assumption of symmetry in cwall
c		
c		Here build up model based on illuminated baffle zones
c
		if(i.ge.88)glaremodel(i,j)=a8*cz*((150.-ty)/150.)		!include Z8 septum
		if(i.ge.86.and.i.lt.88)glaremodel(i,j)=glaremodel(88,j)
		if(i.ge.83.and.i.lt.86)then					!include Z7 septum
			if(ty.lt.80.)then	
			glaremodel(i,j)=glaremodel(86,j)+a7*(float(86-i)/3.)*((95.-ty)/95.)
			else
			glaremodel(i,j)=glaremodel(86,j)
			endif
		endif
		if(i.ge.80.and.i.lt.83)glaremodel(i,j)=glaremodel(83,j)
		if(i.ge.72.and.i.lt.80)then					!include Z6 septum
			if(ty.lt.54.)then
			arg=1.
			if(ty.gt.20.)arg=((110.-ty)/90.)
			glaremodel(i,j)=glaremodel(83,j)+a6*(float(80-i)/8.)*arg
			else
			glaremodel(i,j)=glaremodel(83,j)
			endif
		endif
		if(i.ge.69.and.i.lt.72)glaremodel(i,j)=glaremodel(72,j)
		if(i.ge.67.and.i.lt.69)then					!include Z5 septum
			if(ty.lt.30.)then
			arg=1.
			if(ty.gt.10.)arg=((70.-ty)/60.)
			glaremodel(i,j)=glaremodel(72,j)+a5*(float(69-i)/2.)*arg
			else
			glaremodel(i,j)=glaremodel(72,j)
			endif
		endif
		if(i.ge.63.and.i.lt.67)glaremodel(i,j)=glaremodel(67,j)
		if(i.ge.59.and.i.lt.63)then					!include Z4 septum
			if(ty.lt.25.)then
			arg=1.
			if(ty.gt.8.)arg=((25.-ty)/17.)
			if(i.eq.60.or.i.eq.61)arg=2.*arg
			glaremodel(i,j)=glaremodel(67,j)+a4*arg
			else
			glaremodel(i,j)=glaremodel(67,j)
			endif
		endif
c
c		Here propagate the values below 59 degrees, reducing the septum-bottom contributions appropriately
c		as outer vanes block the view of the bottoms...
c		Note that the Z4 septum bottom is not illuminated below 59 degrees in theta-x
c		 
		if(i.lt.59.and.ty.gt.20.)then
			glaremodel(i,j)=glaremodel(63,j)
			arg1=(58.5-tx)/13.
			if(arg1.gt.1.)arg1=1.
			if(ty.lt.30.)then
			  arg=1.
			  if(ty.gt.10.)arg=((70.-ty)/60.)
			  glaremodel(i,j)=glaremodel(i,j)-a5*arg1*arg		 !take away Z5 septum illumination
			endif
			arg1=(58.5-tx)/70.
			if(ty.lt.54.)then
			  arg=1.
			  if(ty.gt.20.)arg=((110.-ty)/90.)
			  glaremodel(i,j)=glaremodel(i,j)-a6*arg1*arg		 !take away Z6 septum illumination
			endif
			arg1=(58.5-tx)/50.
			if(arg1.gt.1.)arg1=1.
			if(ty.lt.80.)then
			  glaremodel(i,j)=glaremodel(i,j)-a7*arg1*((95.-ty)/95.) !take away Z7 septum illumination
			endif
		endif
c
		glaremodel(i,181-j)=glaremodel(i,j)				 !assumption of symmetry in ty
c		
	enddo
	enddo
	do i=1,90
	do j=1,180
		jj=j-90
		tx=float(i)-0.5
		if(tx.eq.89.5)tx=89.
		ty=float(jj)-0.5
		if(ty.eq.89.5)ty=89.
		if(ty.eq.-89.5)ty=-89.
		cx=sind(tx)/sqrt(1.+((tand(ty))**2) *(cosd(tx)**2))
		cy=sind(ty)/sqrt(1.+((tand(tx))**2) *(cosd(ty)**2))
		arg=cx*cx+cy*cy
		cz=0.
		if(arg.ge.1.0)print *,i,j,arg
		if(arg.lt.1.0)cz=sqrt(1.-arg)
c
c		Here include surface-brightness illumination of septum bottoms
c		Included is a provision for tilting the plane about the Y-axis, to introduce asymmetry
c		
		cx1=s1*cx+c1*cz						!find normal to septum-bottom planes
		cz1=-s1*cz+c1*cx
		cwall=-s2*cy+c2*cx1
c		glaremodel(i,j)=3.*cwall				!handy to view amount of asymmetry, if any
		glaremodel(i,j)=glaremodel(i,j)*cwall
c		glaremodel(i,j)=0.				!**********Kill everything but vane edges
c		
		if(tx.gt.58..or.abs(ty).gt.20.)then
c
c		Here include vane-edge contributions from Z4 through Z7
c		The first value for arg is the length of the illuminated line
c		This is somewhat ad hoc for the long edges since we don't know how far up the edge remains effective
c		The next factor attempts to include how much of the round illuminated edge is in view
c		The third factor, long edges only, is the cosine**2 of the incident x angle for surface brightness
c		This latter: one cos(tx) for amount of light and another for the shadowing of the textured vane edge
c
			arg=1.
			if(abs(ty).gt.10.)arg=(25.-abs(ty))/15.
			arg=arg*(cosd(tx-35.)**2)
			arg=arg*(cosd(tx+25.)**3)	!ad hoc factor to reduce hot spot in Z4 zone near theta-y=0
			if(abs(ty).le.25..and.tx.lt.64.)glaremodel(i,j)=glaremodel(i,j)+z4*arg	!Z4 short edge
			arg=1.
			if(abs(ty).gt.5.)arg=(30.-abs(ty))/25.		
			arg=arg*(cosd(tx-35.)**2)
			arg=arg*cosd(tx-25.)		!ad hoc factor to reduce hot spot in Z5 zone near theta-y=0
			if(abs(ty).le.30..and.tx.lt.69.)glaremodel(i,j)=glaremodel(i,j)+z5*arg	!Z5 short edge
			arg=1.
			if(abs(ty).gt.25.)arg=(54.-abs(ty))/30.		
			arg=arg*(cosd(tx-35.)**2)
			if(abs(ty).lt.54..and.tx.lt.80.)glaremodel(i,j)=glaremodel(i,j)+z6*arg	!Z6 short edge
			arg=1.
			if(abs(ty).gt.10.)arg=(85.-abs(ty))/75.		
			arg=arg*(cosd(tx-35.)**2)
			if(abs(ty).lt.80..and.tx.lt.86.)glaremodel(i,j)=glaremodel(i,j)+z7*arg	!Z7 short edge
c
c		Note that the long edges of Z7 and Z8 are blocked by Z6			
c
			arg=1.
			if(tx.gt.40.)arg=(64.-tx)/24.
			arg=arg*(abs(ty)+40.)/90.
			arg=arg*(cosd(tx-50.)**2)
			arg=arg*(cosd(tx+25)**3)	!ad hoc factor to reduce hot spot in Z4 zone near theta-y=0
			if(abs(ty).le.25..and.tx.lt.64.)glaremodel(i,j)=glaremodel(i,j)+zz4*arg	!Z4 long edge
			arg=1.
			if(tx.gt.50.)arg=(69.-tx)/19.
			arg=arg*(abs(ty)+40.)/90.
			arg=arg*(cosd(tx-35.)**2)
			arg=arg*cosd(tx-25.)		!ad hoc factor to reduce hot spot in Z5 zone near theta-y=0
			if(abs(ty).le.30..and.tx.lt.69.)glaremodel(i,j)=glaremodel(i,j)+zz5*arg	!Z5 long edge
			arg=1.
			if(tx.gt.73.)arg=(80.-tx)/7.
			arg=arg*(abs(ty)+40.)/90.
			arg=arg*(cosd(tx-35.)**2)
			if(abs(ty).lt.54..and.tx.lt.80.)glaremodel(i,j)=glaremodel(i,j)+zz6*arg	!Z6 long edge
c
c		Here include a septum contribution the data call for, that seems to come from some collimation by the 
c		Z3 and Z5 vanes in viewing the Z6 septum: illumination doesn't scatter to the best place behind Z3
c		unless it's pretty far forward towards the rear of Z5....or something like this, ad hoc....
c
c			if(tx.lt.69.and.abs(ty).lt.54.)glaremodel(i,j)=glaremodel(i,j)+b
		endif
c
	enddo
	enddo
	open(22,file='d:\glaremodel\glare_multiplier_model.grd')
	write(22,2)
2	format('DSAA'/'90 180'/'0.5 89.5'/'-89.5 89.5'/'0 3.01')
	do j=1,180
	    write(22,23)(glaremodel(i,j),i=1,90)
23	format(30f5.2)	    
	enddo
	close(22)
	open(22,file='d:\glaremodel\glaredge.grd')
	write(22,2)
	do j=1,180
	    write(22,23)(glaredge(i,j),i=1,90)
	enddo
	close(22)
	stop 'done'
	end