
	Subroutine subtractpolarmaps(ic,instar,histogram,histogramcnt,infile,minorb,totalorbnum,maxmag)
	implicit none
	
	character	dum, starname*15, infile*52, outfile*43, ctype*10
	character	outgrid3*40, outgrid4*37, outgrid7*32
	integer*4	year, day, time, ira, jdec
	integer*4	nxpol, nypol, nxrot, nyrot, direction
	integer*4	b, dloop, dloopmax, bmin, endstar, barg, spec, lum, ic
	integer*4	ibad, k, starcount, imin, imax, jmax, jmin, iteration, xarg, yarg
	integer*4	i, j, ii, jj, xnew, ynew, iarg, jarg, xcen, ycen, xcensave, ycensave
	integer*4	back, zerostar, zerostars, north, median, xcen2, ycen2
	integer*4	badleastchi, badmag
	integer*4	indx(120),m,mzero, kback, kbacksave1,  kbacksave2, zerostarcount
	integer*4	xrot, yrot, npFinal(800,800)
	integer*4	np(800,800), npsub(800,800), npsubcnt(800,800)
	integer*4	instar(200,200), rotstar(200,200), substar(100,100), rotstarcnt(200,200)
	integer*4	newfinal(800,800), npstdall(800,800), npstdallcnt(800,800)
	integer*8	stdstartotal, startotal, chisq, smeitotal
	integer*4	minorb, totalorbnum, histogramcnt(1000,100)
	real*4		tx, ty, px, py, distance, distance0, psi
	real*4		histogram(1000,100), orbitfraction(720,360), fraction, radius2, x2, y2
	real*4		maxmag, raarg, decarg, boxrim(120), rot(720,360), npreal(800,800)
	real*4		secondra, seconddec, secondbrightfactor, chiback, phi0, phi
	real*4		savera, savedec, ra, dec, cra, cdec, x, y, cosq, mag, arg, mymag
	real*4		brightfactor, radius, chi, chimin, tempfactor, smeiadjust, delra, deldec
	real*8		theta, cxpol, cypol, ddpol, cxrot, cyrot, ddrot

	chisq=0
	badmag=0
	badleastchi=0
	
	do i=1,720
	do j=1,360
	  rot(i,j)=0.0
	  orbitfraction(i,j)=0.0
	enddo
	enddo
	
	!!Read in the rotator map & the orbit fraction map
	call smei_sky_read(infile,7,nxrot,nyrot,orbitfraction,cxrot,cyrot,ddrot,ctype)
	call smei_sky_read(infile,8,nxrot,nyrot,rot,cxrot,cyrot,ddrot,ctype)
	
	!!Loop over both polar plots
	do north=1,2		!! 1->north polar plot, 2->south polar plot
	
	do i=1,800
	do j=1,800
	  npstdall(i,j)=0
	  npstdallcnt(i,j)=0
	  np(i,j)=0
	  npFinal(i,j)=0
	  newfinal(i,j)=0
	  npreal(i,j)=0.0
	enddo
	enddo
	
	if(ic.eq.2)then
	  if(north.eq.1)then
	    write(outfile,1245) 'Chi_Output/Chi_Output_c2_',totalorbnum,'_findq_np.txt'
	  else
	    write(outfile,1245) 'Chi_Output/Chi_Output_c2_',totalorbnum,'_findq_sp.txt'
	  endif
	elseif(ic.eq.1)then
	  if(north.eq.1)then
	    write(outfile,1245) 'Chi_Output/Chi_Output_c1_',totalorbnum,'_findq_np.txt'
	  else
	    write(outfile,1245) 'Chi_Output/Chi_Output_c1_',totalorbnum,'_findq_sp.txt'
	  endif
	else
	  if(north.eq.1)then
	    write(outfile,1245) 'Chi_Output/Chi_Output_c3_',totalorbnum,'_findq_np.txt'
	  else
	    write(outfile,1245) 'Chi_Output/Chi_Output_c3_',totalorbnum,'_findq_sp.txt'
	  endif
	endif
1245	format(A25,i5.5,A13)
	open(205,file=outfile,form='formatted')
	write(205,2050)
2050	format('StarName    i2  i1 MyMag  mSMEI  DelMag  RA      Dec     Delra  Del
     &Dec Tx     Ty     phi0 phi  Psi  distance Chi     OrbFraction')
	
	!!Read Polar Plot
	if(north.eq.1)then
	  call smei_sky_read(infile,1,nxpol,nypol,npreal,cxpol,cypol,ddpol,ctype)
	else
	  call smei_sky_read(infile,2,nxpol,nypol,npreal,cxpol,cypol,ddpol,ctype)
	endif
	
	do i=1,800
	do j=1,800
	  np(i,j)=nint(npreal(i,j))
	  npFinal(i,j)=np(i,j)
	enddo
	enddo
	
	open(44,file='BrightStarsSorteda.txt',readonly)
	starcount=0
	zerostarcount=0
	zerostar=0
	zerostars=0
c****** Loop over all bright stars **************************************************
	do k=1,5573
	zerostars=zerostars+zerostar
	zerostar=0
	read(44,444)savera,savedec,mag,spec,lum,starname
444	format(f12.8,2x,f12.8,16x,f5.2,2x,i2,2x,i1,21x,a15)
	if(mag.gt.maxmag) goto 976
	
	!! smeitotal=95.75*4500*(2.5^-mag) this formula must be examined...
c	
c	OK, here's Andy's best take on the #s in this formula: 
c	The 5500 is the number of surface-brightness ADUs that a 0th magnitude star would have assuming 0.55 ADUs is an S10
c	AND (here's the hard part!) the 95.75 is an equivalent factor to the naive 100x one would use at the map's equator,
c	to convert from square degrees to sky bins. 
c	Why not exactly 100? The original scaling from CCD to sky contains radial corrections as in equations (2) and (3) of 
c	my memo dated 11 July 2001, but with the calibrated values and a FOV which is not quite centered on theta-Y = 0 as in
c	the memo. Thus this formula represents the best attempt to get back to ADUs on the CCD, from ADUs on a sky map.
c	On the other hand, we could argue that once sky maps are generated, whether CCD ADUs are perfectly preserved or not,
c	presently they're not, the NOW ADUs are defined to be what we've got on a sky map, which is only approximately that
c	which we see on individual CCD frames. In this case, we should use the 100x. 
c	I'm now (12-10-04) changing this to the	100x. Also, as of this date 0.57 ADUs per S10 looks a bit better...
c
	
	!! Don't do stars that might fall off the edges of np
	if((north.eq.1.and.savedec.lt.40).or.(north.ne.1.and.savedec.gt.-40)) goto 976
	
c	smeitotal=100.*4500.0*(2.5**(-1.0*mag))
	cosq=cosd(savedec)*cosd(savedec)
	direction=0

	do i=1,200
	do j=1,200
	  rotstar(i,j)=0
	  rotstarcnt(i,j)=0
	enddo
	enddo
	
	do i=1,100
	do j=1,100
	  substar(i,j)=0
	enddo
	enddo
	
	!!Find approximate centroid coordinates on sky map from catalogue sky coordinates
	if(north.eq.1)then	
	  xcensave = nint( cxpol+ddpol*(90.0-savedec)*cosd(savera) )
	  ycensave = nint( cypol+ddpol*(90.0-savedec)*sind(savera) )
	else
	  xcensave = nint( cxpol+ddpol*(90.0+savedec)*cosd(savera) )
	  ycensave = nint( cypol+ddpol*(90.0+savedec)*sind(savera) )
	endif
	        
        if(xcensave.lt.25.or.xcensave.gt.775.or.ycensave.lt.25.or.ycensave.gt.775) then
c	  print*, '*** ERROR: Edge of Sky Map ***'
	  ibad=1
	  go to 976
	elseif(np(xcensave,ycensave).eq.0)then
c	  print*, '*** ERROR: No Sky Coverage at Location ***'
	  ibad=1
	  go to 976
	endif
	
	!!Test for holes near edge of map, should never happen now that polar maps are more coarse...
	do i=xcensave-2,xcensave+2
	do j=ycensave-2,ycensave+2
	  if(np(i,j).eq.0)then
c	    print*, '*** ERROR: Holes in Star Location ***'
	    ibad=1
	    go to 976
	  endif
	enddo
	enddo
	
	do i=1,15
	  if(starname(i:i).eq.' '.or.starname(i:i).eq.'*') write(starname(i:i),'(A1)') '_'
	enddo
	
c	print*, '*********************************************************'
c	print 156, k,starname,savera,savedec,mag
c156	format(i4,2x,a16,2x,f7.3,2x,f7.3,2x,f5.2)
c	print*, 'Centroid: ',xcensave, ycensave
	
	imin=xcensave-15
        imax=xcensave+15
        jmin=ycensave-15
        jmax=ycensave+15
	
	  !!Find initial Background using pixels on the rim of the surrounding box
	  m=0
	  mzero=0
	  do i=imin,imax
	  do j=jmin,jmax	    
	    if(i.eq.imin.or.i.eq.imax.or.j.eq.jmin.or.j.eq.jmax)then
	      m=m+1
	      boxrim(m)=float(npFinal(i,j))
	      if(npFinal(i,j).le.0)mzero=mzero+1
	    endif	    
	  enddo
	  enddo
	  
	  !!Sort the pixels on the rim of the surrounding box
	  Call INDEXX(120,boxrim,indx)
	  
	  !!Bail if more than half of the boxrim is outside the map
	  if(mzero.gt.50)then
c	    print*, '*********** Error: Star box outside map ********'
	    ibad=1
	    go to 976
	  endif
	  median=(120+mzero)/2
	  
	  !!Calculate the initial background from the lowest 20 rim pixels
	  back=0
	  do i=mzero+1,mzero+20
	    back=back+int(boxrim(indx(i)))
	  enddo
	  back=back/20
	
	!!Find total surface brightness of star binned in np map
	startotal=0
	do i=imin,imax
	do j=jmin,jmax
	  startotal=startotal+npFinal(i,j)-back	
	enddo
	enddo
	
c	print*, '	 Back(low 20 avg): ',back
c	print*, '	 Startotal: ',startotal
c	print*, '	 Smeitotal[not used]: ',smeitotal
	
	!!Check for bad background if startotal is less than zero
	if(startotal.le.0)then
c	  print*, '	 Lowest Back: ',boxrim(indx(mzero+1)), mzero
	  back=boxrim(indx(mzero+1))
	  startotal=0
	  do i=imin,imax
	  do j=jmin,jmax
	    startotal=startotal+npFinal(i,j)-back	
	  enddo
	  enddo
	  if(startotal.le.0)then
c	    print*, '*********** Error: Bad Background ********'
	    goto 976
	  endif
	endif

	!!Extract the angle of rotation
	xrot = nint( cxrot + ddrot*savera )
	yrot = nint( cyrot + ddrot*savedec )
	theta = rot(xrot,yrot)
	!!Convert the equatorial angle to the proper polar angle
	if(north.eq.1)then
	  theta = 360.0 - (theta+90.0+savera)
	  if(theta.gt.360.0) theta=theta-360.0
	  if(theta.lt.0.0) theta=theta+360.0
	else
	  theta = (theta-savera+90.0)
	  if(theta.gt.360.0) theta=theta-360.0
	  if(theta.lt.0.0) theta=theta+360.0
	endif
	phi0 = theta
	
	!!Find unit vector components in the direction of rotation angle( ~|| velocity vector)
	if(phi0.lt.90.0)then
	  tx = -sind(phi0)
	  ty = -cosd(phi0)
	elseif(phi0.lt.180.0)then
	  ty = sind(phi0-90.0)
	  tx = -cosd(phi0-90.0)
	elseif(phi0.lt.270.0)then
	  tx = sind(phi0-180.0)
	  ty = cosd(phi0-180.0)
	else
	  ty = -sind(phi0-270.0)
	  tx = cosd(phi0-270.0)
	endif
	distance0=sqrt(tx*tx+ty*ty)
	
	!!Stretch and Clockwise Rotate the average star and rebin
	do i=1,200
	do j=1,200	   
	   !!Make the average star finer resolution for rebinning
	   do ii=1,2
	   do jj=1,2	     
	     !!Center of picture is at [0,0]
	     x=float(i)-100.50
	     y=float(j)-100.50
	     
	     if(ii.eq.1) x=x-0.25
	     if(ii.eq.2) x=x+0.25
	     if(jj.eq.1) y=y-0.25
	     if(jj.eq.2) y=y+0.25

	     !!find rotated coordinates
	     xnew=nint( x*dcosd(theta)+y*dsind(theta)+100.5D0 )
	     ynew=nint( -x*dsind(theta)+y*dcosd(theta)+100.5D0 )

	     if(xnew.gt.0.and.xnew.lt.201.and.ynew.gt.0.and.ynew.lt.201)then
	       rotstar(xnew,ynew)=rotstar(xnew,ynew)+instar(i,j)
	       rotstarcnt(xnew,ynew)=rotstarcnt(xnew,ynew)+1
	     endif	     
	   enddo
	   enddo	   
	enddo
	enddo

	!! Average the rotated standard star
	do i=1,200
	do j=1,200
	   if(rotstarcnt(i,j).gt.0) rotstar(i,j)=rotstar(i,j)/rotstarcnt(i,j)
	enddo
	enddo

	!!Put standard star into a smaller box (2.5x2.5 degrees)
	do i=1,200
	do j=1,200
	   xnew=(i-100)+50
	   ynew=(j-100)+50
	   if(xnew.gt.0.and.xnew.lt.101.and.ynew.gt.0.and.ynew.lt.101)then
	     substar(xnew,ynew)=substar(xnew,ynew)+rotstar(i,j)
	   endif
	enddo
	enddo

c*********** Iteration loops to fit least chi value **********************************	
	endstar=0
	iteration=1
	dloopmax=7
c	print*, '  ******* Beginning First Iteration ********'
811	chimin=9999999.0
	if(iteration.eq.2) dloopmax=1		!!don't loop over the direction second time around
	do dloop=1,dloopmax	!! Dec Loop
911	continue
	  
	  do i=1,800
	  do j=1,800
	    npsub(i,j)=0
	    npsubcnt(i,j)=0
	  enddo
	  enddo
	  
	  !!Find pixel coordinates of shifted centroid
	  !!Step size is 0.17 degrees & 10 pixels/degree -> 0.17deg*10pix/deg=1.7pix
	  if(iteration.eq.1)then
	    arg = float(xcensave) + 1.7*float(dloop-4)*tx
	    xcen = nint(arg)
	    arg = float(ycensave) + 1.7*float(dloop-4)*ty
	    ycen = nint(arg)
	  else
	    xcen = xcen2
	    ycen = ycen2
	  endif
	  x=float(xcen)-cxpol
	  y=float(ycen)-cypol
	  
	  !!Find sky coordinates of shifted centroid
	  if(north.eq.1)then	
	    cdec = 90.0 - sqrt(x*x+y*y)/ddpol
	    cra = atan2d(y,x)
	    if(cra.lt.0.0) cra=cra+360.0
	  else
	    cdec = sqrt(x*x+y*y)/ddpol - 90.0
	    cra = atan2d(y,x)
	    if(cra.lt.0.0) cra=cra+360.0
	  endif
	  	  
	  imin=xcen-15
          imax=xcen+15
          jmin=ycen-15
          jmax=ycen+15 
	  
	  !!Put small box star into blank np map
	  do i=1,100
	  do j=1,100
	      x=float(i)-50.50
	      y=float(j)-50.50
	      
	      xarg = xcen + nint(x/4.0)
	      yarg = ycen + nint(y/4.0)

	      if(xarg.gt.0.and.xarg.lt.801.and.yarg.gt.0.and.yarg.lt.801) then
	        npsub(xarg,yarg)=npsub(xarg,yarg)+substar(i,j)
	        npsubcnt(xarg,yarg)=npsubcnt(xarg,yarg)+1
	      endif
	  enddo
	  enddo
	
	  !!Average out the blank map with the avgstar in it, fill boxrim array
	  m=0
	  mzero=0
	  do i=imin,imax
	  do j=jmin,jmax
	    if(npsubcnt(i,j).gt.0) npsub(i,j)=npsub(i,j)/npsubcnt(i,j)	    
	    if(i.eq.imin.or.i.eq.imax.or.j.eq.jmin.or.j.eq.jmax)then
	      m=m+1
	      boxrim(m)=float(npFinal(i,j))
	      if(npFinal(i,j).le.0)mzero=mzero+1
	    endif	    
	  enddo
	  enddo
c	  !!Bail if more than half of the boxrim is outside the map
	  if(mzero.gt.50)then
c	    print*, '*********** Error: Star box outside map ********'
	    ibad=1
	    go to 976
	  endif
	  median=(120+mzero)/2
	  
	  !!Sort the pixels on the rim of the surrounding box
	  Call INDEXX(120,boxrim,indx)
	  
	  !!Background Loop ************************************
	  do kback=0,6
	   if(iteration.eq.1)then
	     back=boxrim(indx(median+(kback-4)*10))
	   else
	     back=boxrim(indx(median+(kbacksave1-4)*10+(kback-3)*2))
	   endif
	   if(endstar.eq.1)then
	     back=boxrim(indx(median+(kbacksave1-4)*10+(kbacksave2-3)*2))
	   endif
	
	   !!Find total surface brightness of the standard star rebinned in np map
	   stdstartotal=0
	   startotal=0
	   do i=imin,imax
	   do j=jmin,jmax
	    if(npsub(i,j).gt.0)then
	      startotal=startotal+npFinal(i,j)-back
	      stdstartotal=stdstartotal+npsub(i,j)
	    endif	
	   enddo
	   enddo
	   if(endstar.eq.1.and.startotal.le.0) then
c	     print*, '****** ERROR: Negative Startotal as Least Chi Value ********'
	     badleastchi=badleastchi+1
	     goto 976
	   elseif(startotal.le.0)then
c	     print*, '****** ERROR: Negative Startotal ********'
	     if(zerostar.eq.0) zerostar=1
	     zerostarcount=zerostarcount+1
	     goto 864
	   endif
	
	   !!Brightness Loop ************************************
	   do b=1,9

	    barg=b
	    if(endstar.eq.1) barg=bmin
	    if(iteration.eq.1) then
	      smeiadjust=1.0+0.1*(barg-5)
	    else
	      smeiadjust=1.0+0.25*(barg-5)
	    endif
	    	    
	    !!Find brightfactor, effectively calibrating the average star
c	    brightfactor=dfloat(smeitotal)/dfloat(stdstartotal)
	    brightfactor=dfloat(startotal)/dfloat(stdstartotal)
	    if(iteration.eq.1)then
	      brightfactor=brightfactor*smeiadjust
	    else
	      brightfactor=secondbrightfactor*smeiadjust
	    endif 
	
	    !!Make final products made only once per star
	    if(endstar.eq.1)then
	    
	      !!Check for crazy magnitude before subtraction, see mymag formula 50 lines below below
	      if(brightfactor.le.0.0) goto 976
	    
	      !!Put the brightfactored standard star into a map to contain all stars, don't average...
	      do i=1,800
	      do j=1,800
	        if(npsub(i,j).gt.0) npstdall(i,j)=npstdall(i,j)+npsub(i,j)*brightfactor
	      enddo
	      enddo
	
	      !!PERFORM SUBTRACTION FROM PRIMARY SKY MAP
	      do i=1,800
	      do j=1,800
	        if(npsub(i,j).gt.0) npFinal(i,j)=nint( float(npFinal(i,j))-float(npsub(i,j))*brightfactor )
	      enddo
	      enddo
	    endif
	
	    !!Calculate the chi value
	    chisq=0
	    if(endstar.eq.0)then
	      do i=imin,imax
	      do j=jmin,jmax
	        arg = float(npFinal(i,j))-float(npsub(i,j))*brightfactor
	        arg = back-arg
	        radius = sqrt( float((i-xcen)*(i-xcen)) + float((j-ycen-2)*(j-ycen-2)) )
	        if(radius.lt.15.)  then
	          chisq=chisq+(arg*arg)
	        endif
	      enddo
	      enddo
	    else
	      do i=imin,imax
	      do j=jmin,jmax
	        arg = npFinal(i,j)
	        arg = back-arg
	        radius = sqrt( float((i-xcen)*(i-xcen)) + float((j-ycen-2)*(j-ycen-2)) )
	        if(radius.lt.15.)  then
	          chisq=chisq+(arg*arg)
	        endif
	      enddo
	      enddo
	    endif
	  
	    chi=sqrt(float(chisq))
	    if(endstar.eq.1)then
	      tempfactor=dfloat(smeitotal)/dfloat(stdstartotal)
	      if(stdstartotal.gt.0.0.and.brightfactor.gt.0.0)then
	        if(north.eq.1)then
	          mymag=-1.0*log(((90.0-savedec)*0.0174533/cosd(savedec))*stdstartotal*brightfactor/450000.0)/log(2.5) !see note above about this formula
	        else
	          mymag=-1.0*log(((90.0+savedec)*0.0174533/cosd(savedec))*stdstartotal*brightfactor/450000.0)/log(2.5) !see note above about this formula	      
	        endif
	      else
	        mymag=-9.9
		badmag=badmag+1
	      endif
	      
cccccccccccccccccccccccccccccc Code from Andy's Starfixum program cccccccccccccccccccccccccccccc

		if(lum.eq.1)mag=mag-0.2		!ad hoc luminosity type adjustment...
		if(lum.eq.2)mag=mag-0.1		!ad hoc luminosity type adjustment...
		if(spec.gt.56)mag=mag-0.25		!ad hoc color adjustments...
		if(spec.ge.51.and.spec.le.56)mag=mag-0.15-float(spec-51)/50.
		if(spec.le.30)mag=mag+0.12    
	      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	      
	      if(north.eq.1)then
	        x = xcensave - 400.5
		y = ycensave - 400.5
		radius = sqrt(x*x+y*y)
		x2 = xcen - 400.5
		y2 = ycen - 400.5
		radius2 = sqrt(x*x+y*y)
		delra = 0.1*radius2*acos((x*x2+y*y2)/sqrt((x*x+y*y)*(x2*x2+y2*y2)))		
	      else
	        x = xcensave - 400.5
		y = ycensave - 400.5
		radius = sqrt(x*x+y*y)
		x2 = xcen - 400.5
		y2 = ycen - 400.5
		radius2 = sqrt(x*x+y*y)
		delra = 0.1*radius2*acos((x*x2+y*y2)/sqrt((x*x+y*y)*(x2*x2+y2*y2)))
	      endif
	      deldec = (radius2 - radius)*0.1
c	      distance=direction*sqrt(delra*delra+deldec*deldec)
c	      distance=direction*sqrt(0.01*(xcen-xcensave)*(xcen-xcensave)+0.01*(ycen-ycensave)*(ycen-ycensave))
	      !!This is the shortest distance between two points on a sphere(great circle segment)
	      distance=direction*acosd(cosd(cdec)*cosd(savedec)*cosd(cra-savera)+sind(cdec)*sind(savedec))
	      
	      !!Find direction angle of the ra/dec shift
	      phi=0.0
	      if(deldec.gt.0.0)then 
	        phi=acosd(abs(deldec)/distance)
	      elseif(deldec.lt.0.0)then
	        if(delra.gt.0.0) phi=asind(abs(delra)/distance) 
	      endif
	      if(cra.gt.0.0.and.phi.ne.0.0)then
	        if(cdec.gt.0.0)then
		  phi=phi
		else
		  phi=phi+90.0
		endif
	      elseif(cra.lt.0.0.and.phi.ne.0.0)then
	        if(cdec.lt.0.0)then
		  phi=phi+180.0
		else
		  phi=phi+270.0
		endif
	      endif
	      
	      !!Extract the fraction of the orbit (use the same coordinates definition as the indexing program)
	      ira = nint( cxrot + ddrot*cra )
	      jdec = nint( cyrot + ddrot*cdec )    
	      fraction=orbitfraction(ira,jdec)
c	      if(fraction.eq.0.0) fraction=9.99
	      
	      !!Find angle between rotation angle[tx,ty] and ra/dec shift angle[delra,deldec]
	      psi = acosd((delra*tx+deldec*ty)/(distance*distance0))  
	      if(distance.lt.0.1) then
	        psi=0.0
	      endif
	      
	      write(205,2051)starname,spec,lum,mymag,mag,mymag-mag,cra,cdec,delra,deldec,
     &	Tx,Ty,nint(phi0),nint(phi),nint(psi),distance,nint(chi),fraction
2051	      format(A10,2x,i2,2x,i1,2x,f5.2,2x,f5.2,2x,f6.2,2x,f6.2,2x,f6.2,2x,f5.2,2x,f5.2,2x,
     &	f5.2,2x,f5.2,2x,i3,2x,i3,2x,i3,2x,f7.3,2x,i6,2x,f7.4)

	      arg = fraction*100.0+0.5
	      jarg = nint(arg)
	      if(jarg.eq.101) jarg=100
	      iarg = totalorbnum-minorb+1
	      if(iarg.gt.0.and.iarg.lt.1001.and.jarg.gt.0.and.jarg.lt.101)then
	        histogram(iarg,jarg)=histogram(iarg,jarg)+distance
	        histogramcnt(iarg,jarg)=histogramcnt(iarg,jarg)+1
	      endif
	      
	      starcount=starcount+1
	      if(mod(starcount,10).eq.0) print*, starcount, ' stars subtracted so far...'
	      goto 976
	    endif
	    if(chi.lt.chimin)then
	      chimin=chi
	      bmin=b
	      kbacksave2=kback
	      if(iteration.eq.1)then
	        xcen2=xcen
	        ycen2=ycen
	        secondbrightfactor=brightfactor
	        kbacksave1=kback
		if(dloop.gt.4) then
		  direction=1
		elseif(dloop.lt.4) then
		  direction=-1
		else
		  direction=0
		endif
	      endif
	    endif
	 
	  enddo !!End b-loop over brightness adjustments
864	  continue
	 enddo  !!End kback-loop over background adjustments
	enddo  !!End d-loop over direction adjustments	  
	if(iteration.eq.1)then
	  iteration=2
c	  print*, '  ******* Beginning Second Iteration ********'
	  goto 811
	endif
	endstar=1
	goto 911
	
976	enddo  !!End k-loop over all bright stars
	close(44)
	
	!!Write blank np Map with oriented standard star in bright star locations
	if(north.eq.1)then
	  write(outgrid3,2001) 'Chi_Output/NpStdStars_c2_',totalorbnum,'_findq.grd'
	else
	  write(outgrid3,2001) 'Chi_Output/SpStdStars_c2_',totalorbnum,'_findq.grd'
	endif
2001	format(A25,i5.5,A10)
	open (22,file = outgrid3, form = 'formatted')
	write(22,201)
201	format('DSAA'/'800 800'/'-40 40'/'-40 40'/'-500 1400')
	do j=1,800
	  write(22,'(20i6)')(npstdall(i,j),i=1,800)
	enddo
	close(22)
	
	!!Write subtracted np Map to a file
	if(north.eq.1)then
	  write(outgrid4,2002) 'Chi_Output/NpFinal_c2_',totalorbnum,'_findq.grd'
	else
	  write(outgrid4,2002) 'Chi_Output/SpFinal_c2_',totalorbnum,'_findq.grd'
	endif
2002	format(A22,i5.5,A10)
	open (23,file = outgrid4, form = 'formatted')
	write(23,201)
	do j=1,800
	  write(23,'(20i6)')(npFinal(i,j),i=1,800)
	enddo
	close(23)
	
	!!Write subtracted np Map to a file
	if(north.eq.1)then
	  write(outgrid7,2005) 'Chi_Output/Np_c2_',totalorbnum,'_findq.grd'
	else
	  write(outgrid7,2005) 'Chi_Output/Sp_c2_',totalorbnum,'_findq.grd'
	endif
2005	format(A17,i5.5,A10)
	open (26,file = outgrid7, form = 'formatted')
	write(26,201)
	do j=1,800
	  write(26,'(20i6)')(np(i,j),i=1,800)
	enddo
	close(26)
	
	close(205)
	
	print*, badmag, ' had crazy magnitudes'
	print*, badleastchi, ' had a negative startotal as least chi value'
	print*, zerostars, ' had negative startotals ',zerostarcount, ' times'
	print*, starcount, 'Stars Subtracted'
	
	enddo !!End north loop over both polar plots

	End

c** This function is defined in the analagous SubtractEqMaps.f file
c*********************************************************************************
c      SUBROUTINE INDEXX(N,ARRIN,INDX)	!from Numerical Receipes...
c      implicit none
c      integer*4 N
c      integer*4 INDX(N),I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
c      real*4 ARRIN(N),Q
c      DO 11 J=1,N
c        INDX(J)=J
c11    CONTINUE
c      L=N/2+1
c      IR=N
c10    CONTINUE
c        IF(L.GT.1)THEN
c          L=L-1
c          INDXT=INDX(L)
c          Q=ARRIN(INDXT)
c        ELSE
c          INDXT=INDX(IR)
c          Q=ARRIN(INDXT)
c          INDX(IR)=INDX(1)
c          IR=IR-1
c          IF(IR.EQ.1)THEN
c            INDX(1)=INDXT
c            RETURN
c          ENDIF
c        ENDIF
c        I=L
c        J=L+L
c20      IF(J.LE.IR)THEN
c          IF(J.LT.IR)THEN
c            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
c          ENDIF
c          IF(Q.LT.ARRIN(INDX(J)))THEN
c            INDX(I)=INDX(J)
c            I=J
c            J=J+J
c          ELSE
c            J=IR+1
c          ENDIF
c        GO TO 20
c        ENDIF
c        INDX(I)=INDXT
c      GO TO 10
c      END
c*********************************************************************************

