C+
C NAME:
C	avgstar_subtract_brx
C MODIFICATION HISTORY:
C	July, 2003 A. Smith (UCSD/CASS)
C-
	program avgstar_subtract_brx
	implicit none

	character	strfile*57, stdfile*57, infile3*30, stroutfile*45, stdoutfile*45, scatfile*30
	character	strlist*26 , stdlist*26, dum*20, starname*30
	integer*4	str(200,200), std(200,200), avgstar(200,200)
	integer*4       i, j,ss, date, orbitnum, maxr, minr, x0, y01, y02, ii, jj, back
	integer*4	imin, imax, jmin, jmax, iarg, jarg, spec, lum, avgbackcnt, numsubpixels
	real*4		strres, stdres, mag1, mag2, mag3, elipser1, elipser2
	real*8		strtotal, stdtotal, avgtotal, sqztotal, thetax, strmult, stdmult, smeimult
	real*8		totalres, totalback, resfraction
	real*8		radius, mult, substar(200,200), sqzsubstar(200,200), sqzavgstar(200,200)
	real*8		total, tempavgstar(800,800), ftempavgstar(800,800), rad1, rad2, multdiscrep

	print*, 'Fortran in...'

	!!Read in average star to be subtracted from star transits
	  write (infile3,'(A30)') 'Star_Transits/avg_star_str.grd'
	  open(11,file=infile3,form='formatted',readonly)
	  read(11,119)dum
	  read(11,119)dum
	  read(11,119)dum
	  read(11,119)dum
	  read(11,119)dum
119	  format(a20)
	  read(11,11)((avgstar(i,j),i=1,200),j=1,200)
11	  format(20i6)
	  close(11)
	
	write(strlist,'(A26)') 'Star_Transits/str_list.txt'
	open(18,file=strlist,status='old')
	write(stdlist,'(A26)') 'Star_Transits/std_list.txt'
	open(19,file=stdlist,status='old')
	write(scatfile,'(A30)') 'Star_Transits/res_fraction.dat'
	open(22,file=scatfile,form='formatted')

	y01=106
	elipser1=1.6
	rad1=37.D0
	y02=104
	elipser2=1.75
	rad2=35.D0
	multdiscrep=0.D0
	totalres=0.D0

	do ss=1,27
	  read(18,18) starname
18	  format(A30)
	  write(strfile,'(A27,A30)') 'Star_Transits/str_transits/',starname
	  read(19,19) starname, thetax, mag1, mag2, mag3, spec, lum
19	  format(A30,2x,F10.6,2x,F5.2,2x,F5.2,2x,F5.2,2x,I2.2,2x,I2.2)
	  if(dabs(thetax).lt.10.D0)then
	    rad2=37
	  elseif(dabs(thetax).lt.20.D0)then
	    rad2=36
	  else
	    rad2=35
	  endif
	  !! ad hoc color/brightness adjustments for subtraction
	  if(lum.eq.1)then
	    mag3=mag3-0.15
	  elseif(lum.eq.2)then
	    mag3=mag3-0.15
	  elseif(lum.eq.3)then
	    mag3=mag3-0.1
	  elseif(lum.eq.4)then
	    mag3=mag3-0.1
	  elseif(lum.eq.5)then
	    mag3=mag3-0.05
	  endif
	  if(spec.gt.56)then
	    mag3=mag3-0.40
	  elseif(spec.le.30)then
	    mag3=mag3+0.20
	  else
	    mag3=mag3-0.20-(float(spec-51)/50.D0)
	  endif


	  if(ss.gt.0)then

	  write(stdfile,'(A27,A30)') 'Star_Transits/std_transits/',starname
	  write(stroutfile,'(A18,A15,A12)') 'Star_Transits/sub/',starname(1:15),'_str_sub.grd'
	  write(stdoutfile,'(A18,A15,A12)') 'Star_Transits/sub/',starname(1:15),'_std_sub.grd'
	  
	  do i=1,800
	    do j=1,800
	      tempavgstar(i,j)=0.D0
	      ftempavgstar(i,j)=0.D0
	    enddo
	  enddo

	  do i=1,200
	    do j=1,200
	      str(i,j)=0
	      std(i,j)=0
	      substar(i,j)=0.D0
	      sqzsubstar(i,j)=0.D0
	      sqzavgstar(i,j)=0.D0
	    enddo
	  enddo
	
	!! put avgstar into much finer grid, maintaining total response
	  do i=1,200
	    do j=1,200
	      imax=i*4
	      jmax=j*4
	      imin=imax-3
	      jmin=jmax-3
	      do ii=imin,imax
	        do jj=jmin,jmax
	          tempavgstar(ii,jj)=float(avgstar(i,j))/16.D0
	        enddo
	      enddo
	    enddo
	  enddo

	!! squeeze the fine grid avg star
	  do i=1,800
	    do j=1,800
	      ii = nint( (float(i)-400.5)*dcosd(thetax) + 400.5 )
	      jj = j
	      ftempavgstar(ii,jj)=tempavgstar(i,j)
	    enddo
	  enddo

	!! rebin the fine grid to the normal size grid
	  do i=1,200
	    do j=1,200
	      imax=i*4
	      jmax=j*4
	      imin=imax-3
	      jmin=jmax-3
	      do ii=imin,imax
	        do jj=jmin,jmax
	          sqzavgstar(i,j)=sqzavgstar(i,j)+ftempavgstar(ii,jj)
	        enddo
	      enddo
	    enddo
	  enddo

	  !! read in headers, then image array
	  open(12,file=strfile,form='formatted',readonly)
	  read(12,120)dum
	  read(12,120)dum
	  read(12,120)dum
	  read(12,120)dum
	  read(12,120)dum
120	  format(a20)
	  read(12,12)((str(i,j),i=1,200),j=1,200)
12	  format(20i6)
	  close(12)
	  
	  !! read in headers, then image array
	  open(13,file=stdfile,form='formatted',readonly)
	  read(13,130)dum
	  read(13,130)dum
	  read(13,130)dum
	  read(13,130)dum
	  read(13,130)dum
130	  format(a20)
	  read(13,13)((std(i,j),i=1,200),j=1,200)
13	  format(20i6)
	  close(13)
	  
	  back=65000
	  numsubpixels=0
	  do i=1,200
	    do j=1,200
	      radius = sqrt( float((i-100)*(i-100))*elipser2 + float((j-y02)*(j-y02)) )
	      if(radius.lt.rad2.and.std(i,j).gt.0.and.std(i,j).lt.back) back=std(i,j)
	    enddo
	  enddo

	  strtotal=0.D0
	  avgtotal=0.D0
	  stdtotal=0.D0
	  sqztotal=0.D0
	  do i=1,200
	    do j=1,200
	      radius = sqrt( float((i-100)*(i-100))*elipser1 + float((j-y01)*(j-y01)) )
	      if(radius.lt.rad1.and.str(i,j).gt.0)then
		  strtotal=strtotal+float(str(i,j))
		  avgtotal=avgtotal+float(avgstar(i,j))
	      endif
	      
	      radius = sqrt( float((i-100)*(i-100))*elipser2 + float((j-y02)*(j-y02)) )
	      if(radius.lt.rad2.and.std(i,j).gt.0)then
		  stdtotal=stdtotal+float(std(i,j))
		  sqztotal=sqztotal+sqzavgstar(i,j)
		  numsubpixels=numsubpixels+1
	      endif
	    enddo
	  enddo
	  totalback=float(numsubpixels*back)

	  strmult=strtotal/avgtotal
	  stdmult=stdtotal/sqztotal
	  smeimult=0.2*(2.55**(-mag3))
	  print*, starname

	  sqztotal=0.D0
	  do i=1,200
	    do j=1,200
	      radius = sqrt( float((i-100)*(i-100))*elipser1 + float((j-y01)*(j-y01)) )
	      if(radius.lt.rad1)then
	        if(str(i,j).gt.0)  substar(i,j)=float(str(i,j)) - float(avgstar(i,j))*strmult
	      else
	        substar(i,j)=float(str(i,j))
	      endif

	      radius = sqrt( float((i-100)*(i-100))*elipser2 + float((j-y02)*(j-y02)) )
	      if(radius.lt.rad2)then
	        if(std(i,j).gt.0) then
		  sqzsubstar(i,j)=float(std(i,j)) - sqzavgstar(i,j)*smeimult
		  sqztotal=sqztotal+(sqzavgstar(i,j)*smeimult)
		endif
	       else
		sqzsubstar(i,j)=float(std(i,j))
	      endif
	    enddo
	  enddo

	  strres=0.0
	  stdres=0.0
	  do i=1,200
	    do j=1,200
	      radius = sqrt( float((i-100)*(i-100))*elipser1 + float((j-y01)*(j-y01)) )
	      if(radius.lt.rad1)then
	        strres=strres+substar(i,j)
	      endif

	      radius = sqrt( float((i-100)*(i-100))*elipser2 + float((j-y02)*(j-y02)) )
	      if(radius.lt.rad2.and.std(i,j).gt.0)then
	        stdres=stdres+sqzsubstar(i,j)
	      endif
	    enddo
	  enddo
C	  print*, 'Residues: ',strres, stdres
	  stdres=stdres-totalback
	  resfraction=stdres/sqztotal
	  print*, '			',resfraction

	  !! write out data file, not including double stars
	  if(ss.ne.3.and.ss.ne.4.and.ss.ne.25.and.ss.ne.26)then
	    write(22,22) mag1, mag2, mag3, spec, lum, resfraction, starname(1:14)
22	    format(F5.2,2x,F5.2,2x,F5.2,2x,I2.2,2x,I2.2,2x,F10.6,2x,A14)
	  endif

	  !!Write subtracted image to file
	  open (10,file = stroutfile, form = 'formatted')
	  write(10,10)
10	  format('DSAA'/'200 200'/'0 200'/'0 200'/'-5000 5000')
	  do j=1,200
	    write(10,'(20i6)')(nint(substar(i,j)),i=1,200)
	  enddo
	  close(10)

	  !!Write subtracted image to file
	  open (13,file = stdoutfile, form = 'formatted')
	  write(13,23)
23	  format('DSAA'/'200 200'/'0 200'/'0 200'/'-5000 5000')
	  do j=1,200
	    write(13,'(20i6)')(nint(sqzsubstar(i,j)),i=1,200)
	  enddo
	  close(13)

	  endif
	enddo !!end ss do loop to include all the stars
	close(18)
	close(19)
	close(22)

	print*, 'Fortran out!!!'

	end
