
	Subroutine subtracteqmaps(ic,instar,histogram,histogramcnt,infile,minorb,totalorbnum,maxmag)
	implicit none
	
	character	dum, starname*15, infile*52, outfile1*43, ctype*10
	character	outgrid3*40, outgrid4*37, outgrid7*32
	integer*4	nxeq, nyeq, nxrot, nyrot
	integer*4	year, day, time, dloopmax, ira, jdec, direction
	integer*4	b, dloop, bmin, endstar, barg, spec, lum, ic
	integer*4	k, starcount, imin, imax, jmax, jmin, iteration, xarg, yarg
	integer*4	i, j, ii, jj, xnew, ynew, iarg, jarg, xcen, ycen, back
	integer*4	zerostar, zerostars, median, badmag, badleastchi
	integer*4	indx(120),m,mzero, kback, kbacksave1,  kbacksave2, zerostarcount
	integer*4	eq2Final(3600,1200)
	integer*4	eq2(3600,1200), eqsub(3600,1200), eqsubcnt(3600,1200)
	integer*4	instar(200,200), rotstar(200,200), substar(100,100), rotstarcnt(200,200)
	integer*4	eqstdall(3600,1200), eqstdallcnt(3600,1200)
	integer*4	minorb, totalorbnum, histogramcnt(1000,100)
	integer*8	stdstartotal, startotal, chisq, smeitotal
	real*4		phi0, phi, tx, ty, px, py, distance, distance0, psi
	real*4		eqreal(3600,1200), rot2(720,360), histogram(1000,100), orbitfraction(760,320)
	real*4		boxrim(120), secondra, seconddec, secondbrightfactor, chiback, maxmag, fraction
	real*4		savera, savedec, ra, dec, cra, cdec, x, y, cosq, mag, arg, mymag, fxmag
	real*4		brightfactor, radius, chi, chimin, tempfactor, smeiadjust, delra, deldec
	real*8		theta, cxeq, cyeq, ddeq, cxrot, cyrot, ddrot
	
	chisq=0
	badmag=0
	badleastchi=0
	
	do i=1,720
	do j=1,360
	  rot2(i,j)=0.0
	  orbitfraction(i,j)=0.0
	enddo
	enddo
	
	do i=1,3600
	do j=1,1200
	  eqstdall(i,j)=0
	  eqstdallcnt(i,j)=0
	  eq2(i,j)=0
	  eq2Final(i,j)=0
	enddo
	enddo
	
1234	format(A25,i5.5,A13)
	if(ic.eq.1)then
	  write(outfile1,1234) 'Chi_Output/Chi_Output_c1_',totalorbnum,'_findq_eq.txt'
	elseif(ic.eq.2)then
	  write(outfile1,1234) 'Chi_Output/Chi_Output_c2_',totalorbnum,'_findq_eq.txt'
	else
	  write(outfile1,1234) 'Chi_Output/Chi_Output_c3_',totalorbnum,'_findq_eq.txt'
	endif
	open(105,file = outfile1 ,form='formatted')
	write(105,1050)
1050	format('StarName    i2  i1 MyMag  mSMEI  DelMag  RA      Dec     Delra  Del
     &Dec Tx     Ty     phi0 phi  Psi  distance Chi     OrbFraction')

	call smei_sky_read(infile,0,nxeq,nyeq,eqreal,cxeq,cyeq,ddeq,ctype)
	call smei_sky_read(infile,7,nxrot,nyrot,orbitfraction,cxrot,cyrot,ddrot,ctype)
	call smei_sky_read(infile,8,nxrot,nyrot,rot2,cxrot,cyrot,ddrot,ctype)
	
	do i=1,3600
	do j=1,1200
	   eq2(i,j)=nint(eqreal(i,j))
	   eq2Final(i,j)=nint(eqreal(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*5500*(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	smeitotal=100.*4500.0*(2.5**(-1.0*mag))
c	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
	xcen = nint( cxeq + ddeq*savera )
	ycen = nint( cyeq + ddeq*savedec )        
        if(xcen.lt.25.or.xcen.gt.3575.or.ycen.lt.25.or.ycen.gt.1175) then
c	  print*, '*** ERROR: Edge of Sky Map ***'
	  go to 976
	elseif(eq2(xcen,ycen).eq.0)then
c	  print*, '*** ERROR: No Sky Coverage at Location ***'
	  go to 976
	endif
	
	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)
	
        imin=xcen-15
        imax=xcen+15
        jmin=ycen-15
        jmax=ycen+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(eq2Final(i,j))
	      if(eq2Final(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 ********'
	    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 Eq map
	startotal=0
	do i=imin,imax
	do j=jmin,jmax
	  startotal=startotal+eq2Final(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.le.0)then
	  startotal=0
	  back=boxrim(indx(mzero+1))
c	  print*, '	 Lowest Back: ',back, mzero
	  do i=imin,imax
	  do j=jmin,jmax
	    startotal=startotal+eq2Final(i,j)-back	
	  enddo
	  enddo
	  if(startotal.le.0)then
c	    print*, '*********** Error: Bad Background ********'
	    goto 976
	  endif
	endif

	!!Extract the angle of rotation, Paul's maps are righthanded(unlike Nortons) -> reflect angle about y-axis
	xarg = nint( cxrot + ddrot*savera )
	yarg = nint( cyrot + ddrot*savedec ) 
	theta = 360.0 - rot2(xarg,yarg)
	if(theta.gt.360.0) theta = theta - 360.0
	if(theta.lt.0.0) theta = theta + 360.0
	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 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	!! Direction Loop
911	continue
	  
	  do i=1,3600
	  do j=1,1200
	    eqsub(i,j)=0
	    eqsubcnt(i,j)=0
	  enddo
	  enddo	  
	  
	  if(iteration.eq.1)then
	    cra=savera+0.17*float(dloop-4)*tx
	    cdec=savedec+0.17*float(dloop-4)*ty
	  else
	    cra=secondra
	    cdec=seconddec
	  endif
	  arg=cosd(cdec)
	  cosq=arg*arg
	  
	  !!Find shifted centroid coordinates
	  xcen = nint( cxeq + ddeq*savera )
	  ycen = nint( cyeq + ddeq*savedec )
	  
	  imin=xcen-15
          imax=xcen+15
          jmin=ycen-15
          jmax=ycen+15 
	  
	  !!Put small box star into blank eq map
	  do i=1,100
	  do j=1,100
	    x=float(i)-50.50
	    y=float(j)-50.50
	  
	    radius = sqrt(x*x*cosq + (y-3.0)*(y-3.0))
	  
	    if(radius.lt.40.0)then
	      dec=y*(2.5/100.0)+cdec
	      ra=x*(2.5/100.0)/cosd(dec)+cra

	      xarg = nint( cxeq + ddeq*ra )
	      yarg = nint( cyeq + ddeq*dec )

	      if(xarg.gt.0.and.xarg.lt.3601.and.yarg.gt.0.and.yarg.lt.1201) then
	        eqsub(xarg,yarg)=eqsub(xarg,yarg)+substar(i,j)
	        eqsubcnt(xarg,yarg)=eqsubcnt(xarg,yarg)+1
	      endif
	    endif
	  enddo
	  enddo
	  
	  !!Average out the blank map with the avgstar in it
	  do i=imin,imax
	  do j=jmin,jmax
	    if(eqsubcnt(i,j).gt.0) eqsub(i,j)=eqsub(i,j)/eqsubcnt(i,j)	    	    
	  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(eqsubcnt(i,j).gt.0) eqsub(i,j)=eqsub(i,j)/eqsubcnt(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(eq2Final(i,j))
	      if(eq2Final(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 ********'
	    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
c	     back=boxrim(indx(mzero+1+10*kback))
	     back=boxrim(indx(median+(kback-4)*10))
	   else
c	     back=boxrim(indx(mzero+1+10*kbacksave1+2*(kback-3)))
	     back=boxrim(indx(median+(kbacksave1-4)*10+(kback-3)*2))
	   endif
	   if(endstar.eq.1)then
c	     back=boxrim(indx(mzero+1+10*kbacksave1+2*(kbacksave2-3)))
	     back=boxrim(indx(median+(kbacksave1-4)*10+(kbacksave2-3)*2))
	   endif
	
	   !!Find total surface brightness of the standard star rebinned in Eq map
	   stdstartotal=0
	   startotal=0
	   do i=imin,imax
	   do j=jmin,jmax
	    if(eqsub(i,j).gt.0)then
	      startotal=startotal+eq2Final(i,j)-back
	      stdstartotal=stdstartotal+eqsub(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
	     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,3600
	      do j=1,1200
	        if(eqsub(i,j).gt.0) eqstdall(i,j)=eqstdall(i,j)+eqsub(i,j)*brightfactor
	      enddo
	      enddo
	
	      !!PERFORM SUBTRACTION FROM PRIMARY SKY MAP
	      do i=1,3600
	      do j=1,1200
	        if(eqsub(i,j).gt.0) eq2Final(i,j)=nint( float(eq2Final(i,j))-float(eqsub(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(eq2Final(i,j))-float(eqsub(i,j))*brightfactor
	        arg = back-arg
	        radius = sqrt( float((i-xcen)*(i-xcen))*cosq + float((j-ycen-2)*(j-ycen-2)) )
	        if(radius.lt.14)  then
	          chisq=chisq+(arg*arg)
	        endif
	      enddo
	      enddo
	    else
	      do i=imin,imax
	      do j=jmin,jmax
	        arg = eq2Final(i,j)
	        arg = back-arg
	        radius = sqrt( float((i-xcen)*(i-xcen))*cosq + float((j-ycen-2)*(j-ycen-2)) )
	        if(radius.lt.14)  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
	        mymag=-1.0*log(cosd(cdec)*stdstartotal*brightfactor/450000.0)/log(2.5)	!see note above about this formula
	      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	       
	      
	      delra = cosd(cdec)*(cra-savera)
	      deldec = cdec-savedec
	      distance=direction*sqrt(delra*delra+deldec*deldec)
	      
	      !!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
	      xarg = nint( cxrot + ddrot*cra )
	      yarg = nint( cyrot + ddrot*cdec )    
	      fraction=orbitfraction(xarg,yarg)
c	      if(fraction.eq.0.0) fraction=9.99	!!Not sure why this is here
	      
	      !!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(105,1051)starname,spec,lum,mymag,mag,mymag-mag,cra,cdec,delra,deldec,
     &	Tx,Ty,nint(phi0),nint(phi),nint(psi),distance,nint(chi),fraction
1051	      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 checked so far...'
	      goto 976
	    endif
	    if(chi.lt.chimin)then
	      chimin=chi
	      bmin=b
	      kbacksave2=kback
	      if(iteration.eq.1)then
	        secondra=cra
	        seconddec=cdec
	        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 adjustment	  
	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 EQ Map with oriented standard star in bright star locations
2001	format(A25,i5.5,A10)
	write(outgrid3,2001) 'Chi_Output/EqStdStars_c2_',totalorbnum,'_findq.grd'
	open (22,file = outgrid3, form = 'formatted')
	write(22,202)
	do j=1,1200
	  write(22,'(20i6)')(eqstdall(i,j),i=1,3600)
	enddo
	close(22)
	
	!!Write subtracted EQ Map to a file
2002	format(A22,i5.5,A10)
	write(outgrid4,2002) 'Chi_Output/EqFinal_c2_',totalorbnum,'_findq.grd'
	open (23,file = outgrid4, form = 'formatted')
	write(23,202)
202	format('DSAA'/'3600 1200'/'0 360'/'-60 60'/'-500 1400')
	do j=1,1200
	  write(23,'(20i6)')(eq2Final(i,j),i=1,3600)
	enddo
	close(23)
	
	!!Write EQ Map to a file
2005	format(A17,i5.5,A10)
	write(outgrid7,2005) 'Chi_Output/Eq_c2_',totalorbnum,'_findq.grd'
	open (26,file = outgrid7, form = 'formatted')
	write(26,202)
	do j=1,1200
	  write(26,'(20i6)')(eq2(i,j),i=1,3600)
	enddo
	close(26)
	
	close(105)
	
	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'

	End

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

