c+
c name:		star2smei.for
c purpose:	
c	This program reads SIMBAD Query Result files, and decodes the Spectral & Luminosity Class
c	A "SMEI magnitude" is calculated from Mv and a correction table generated by program "MV2SMEI.for"
c	Objects are listed, but not included in the output file, which have no spectral type,
c	unless they are clusters.
c
c	Spectral types "C,S,N,R" are assigned respectively K5,M4,M4,K0, all Lum III
c	Wolf-Rayet stars are assigned K2 or M3 for respective spectral class "C" or "M", all Lum III
c	Subtypes, if unavailable, are assigned = 5
c	Luminosities, if unavailable, are assigned to III (G & redder) or V (otherwise)
c
c inputs:	stars2smei_filled.txt and 50 SIMBAD files covering the whole sky to 8th vis mag
c outputs:	
c	The output file is created which begins with the spectral type field, followed by otype, 
c	then RA, Dec, the three magnitudes, the found or assigned spectral class,
c	and finally the identifier field. This new file requires that Msmei be brighter than 6th.
c calls:	subroutine 'double(ifirst,identifier,iout)'
c mods:		Jan 2003, Andrew Buffington (UCSD; (858)-534-6630; abuffington@ucsd.edu)
c-
	program star2smei
	implicit none
c
	integer*4 i,j,k,irahour,iramin,idecdeg,idecmin,idecsec,iarg,jarg,icluster,icluscnt,iwolf
	integer*4 icount,ibad,igood,namesv,numbsv,icnt,ilumcnt,jnumbsv,ikeep,is,iw,iout,ioutcnt
	integer*4 ilum,lumok,ibump,lumsv,itype(70,5),iflag,iflagcnt,iibad,iscnt,m,kmax
	real*4 amagb,amagv,subfac,arg,rasec,asubtype,mv2smei(70,5),amagsmei
	character*130 heads(22)
	character*12 starsin
	character*13 starsout
	character*21 filledin
	character*26 identifier
	character*4 otype,cluster(3),eclipse(4)
	character*1 obafg,subtype,lum(6),names(13),numbs(10),lums(10)
	data cluster /' Cl*',' OpC',' GlC'/
	data eclipse /' EB*',' Al*',' bL*',' WU*'/
	data names /'O','B','A','F','G','K','M',' ','C','S','W','N','R'/
	data numbs /'0','1','2','3','4','5','6','7','8','9'/	!'m','.',' ' are also seen
	data lums  /'I','V','.',' ','m','s','p','e',' ',':'/
	data itype /350*0/
c	
	filledin='stars2smei_filled.txt'
	open (10,file=filledin,readonly)
	do i=1,70
	read(10,10)(mv2smei(i,j),j=1,5)
10	format(5x,5f8.3)
	enddo
	close(10)
	call double(1,identifier,iout)
	icount=0
	iflagcnt=0
	ilumcnt=0
	iscnt=0
	iibad=0
	icluscnt=0
	iwolf=0
	ikeep=0
	is=0
	iw=0
	ioutcnt=0
	do m=1,4
	   kmax=0
	   if(m.eq.2.or.m.eq.3)kmax=23
c
	   if(m.eq.1)starsin = 'n__pole.txt'
	   if(m.eq.2)starsin = 'field+xx.txt'
	   if(m.eq.3)starsin = 'field-xx.txt'
	   if(m.eq.4)starsin = 's__pole.txt'
	   if(m.eq.1)starsout = 'n_poleout.dat'
	   if(m.eq.2)starsout = 'fieldout+.dat' 
	   if(m.eq.3)starsout = 'fieldout-.dat' 
	   if(m.eq.4)starsout = 's_poleout.dat'
	   open (11,file=starsout)
	   do k=0,kmax
		if(m.eq.2.or.m.eq.3)write(starsin(7:8),'(i2.2)')k
		open (12,file=starsin,readonly)
		read (12,14)heads
14		format(a130)
		icnt=0
		do i=1,10000
		 read (12,13)identifier,otype,irahour,iramin,rasec,idecdeg,idecmin,idecsec,amagb,amagv,obafg,subtype,lum
13		 format(a26,1x,a4,2x,i2,1x,i2,1x,f5.2,2x,i3,1x,i2,1x,i2,8x,f5.2,6x,f5.2,3x,8a1,20x)
			if(irahour.eq.0.and.iramin.eq.0.and.rasec.eq.0.0)go to 99	!detects EOF
			if(amagb.eq.0.00000)go to 98			!delete "IRAS variables" &c
			if(m.eq.1.and.idecdeg.lt.60)go to 98					!N pole
			if(m.eq.2.and.(irahour.ne.k.or.idecdeg.lt.0.or.idecdeg.ge.60))go to 98	!N<60deg
			if(m.eq.3.and.(irahour.ne.k.or.idecdeg.gt.0.or.idecdeg.le.-60))go to 98	!S>-60deg
			if(m.eq.4.and.idecdeg.gt.-60)go to 98					!S pole
c			
			arg=2.5**(-amagv)
			subfac=2.5*arg
			icount=icount+1
			icnt=icnt+1
			ibad=1
			igood=0
			lumok=0
			ilum=0
			lumsv=0
			namesv=0
			numbsv=0
			ibump=0
			iflag=0
			icluster=0
			do j=1,10
			  if(j.le.7.and.obafg.eq.names(j))then
				ibad=0
				namesv=j
			  endif
			  if(subtype.eq.numbs(j))then
				igood=1
				numbsv=j
				asubtype=float(j-1)
			  endif
			  if(j.le.8.and.lum(1).eq.lums(j))then
				lumok=1
				lumsv=j
			  endif
			enddo
c			
			if(obafg.eq.names(9))then	!For "C": use K5III 
			  namesv=6
			  asubtype=5.
			  ilum=3
			  iflag=3
			  go to 97
			endif
			if(obafg.eq.names(10))then	!For "S": use M4III 
			  namesv=7
			  asubtype=4.
			  ilum=3
			  iflag=3
			  go to 97
			endif
			if(obafg.eq.names(11))then	!For Wolf-Rayet stars
			   if(subtype.eq.names(9).or.subtype.eq.names(12))then
			   	    ilum=3
			   	    iflag=4
			   	    if(subtype.eq.names(9))then	!C: call it K2III
			   		namesv=6
			   		asubtype=2.
			   	    else			!N: call it M3III
			   	    	namesv=7
			   	    	asubtype=3.
			   	    endif
			   	    go to 97
			   endif
			endif
			if(obafg.eq.names(12))then	!For "N": use M4III 
			  namesv=7
			  asubtype=4.
			  ilum=3
			  iflag=3
			  go to 97
			endif
			if(obafg.eq.names(13))then	!For "R": use K0III 
			  namesv=6
			  asubtype=0.
			  ilum=3
			  iflag=3
			  go to 97
			endif
c
c	Special case for subtype 10: we haven't seen one yet, but maybe...???
c			
			if(subtype.eq.numbs(2).and.lum(1).eq.numbs(1))then
			  print 1000,icnt,otype,obafg,subtype,lum
1000			  format(i6,2x,a4,2x,', Subtype 10!!!    ',8a1)
			  go to 98
			endif
c
c	Clusters...use baseline B5 for open or general clusters, use F8 for globulars
c
			if(ibad.eq.1)then
			   if(otype.eq.cluster(1).or.otype.eq.cluster(2))then
			   	ibad=0
			   	igood=1
			   	icluster=1
			   	namesv=2
			   	asubtype=5.
			   	ilum=5
			   endif
			   if(otype.eq.cluster(3))then
			   	ibad=0
			   	igood=1
			   	icluster=1
			   	namesv=4
			   	asubtype=5.
			   	ilum=5
			   endif
			endif
c
c			Accept only 7 spectral types: otherwise bail...
c			
			if(ibad.eq.1)then
c			   print 1001,icnt,identifier,otype,obafg,subtype,lum,amagv
1001			   format(i6,2x,a26,a4,2x,', Bad Name    =    ',8a1,' Mv = ',f5.2)
			   iibad=iibad+1
			   go to 98
			endif
c			
			if(igood.ne.1)then		!if subtype not 0->9, just assign it & luminosity
			   asubtype=5.
			   if(namesv.le.4)ilum=5
			   if(namesv.ge.5)ilum=3
			   iflag=1
c			   print 1002,icnt,identifier,otype,obafg,subtype,lum
1002			   format(i6,2x,a26,a4,2x,', Subtype  -> 5 ',3x,8a1)
			   go to 97
			endif			
c
c			Only arrive here if both spectral type & subtype OK
c			
			  if(lumsv.eq.3)then			!check for "#.#", add 0.# if found
			     jnumbsv=0
			     do j=2,9
			     	if(lum(2).eq.numbs(j))jnumbsv=j
			     enddo
			     if(jnumbsv.ne.0)then
			        asubtype=asubtype+0.1*float(jnumbsv)
			        ibump=2
			     else
			     	if(lum(3).eq.lums(3))then	!assign luminosity when "..."
			   		if(namesv.le.4)ilum=5
			   		if(namesv.ge.5)ilum=3
			   		go to 97
			   	else				!bail when neither #.# nor "..."
			     		print 1003,icnt,identifier,otype,obafg,subtype,lum
1003			     		format(i6,2x,a26,a4,2x,', Bad Lums    =    ',8a1)
			     		iibad=iibad+1
			     		go to 98
			     	endif
			     endif
			  endif
c			  
			  if(lum(1+ibump).eq.lums(10))then  !find colon after subtype #, skip it
			     lumok=1
c			     print 1004,icnt,identifier,otype,obafg,subtype,lum
1004			     format(i6,2x,a26,a4,2x,', Find colon...    ',8a1)
			     ibump=ibump+1
			  endif
c
			  if(lumok.ne.1.or.lumsv.gt.3)then		!any other trouble, call it III or V
			    iflag=2
			    if(namesv.le.4)then
			    ilum=5
c			    print 1005,icnt,identifier,otype,obafg,subtype,lum
1005			    format(i6,2x,a26,a4,2x,', Luminosity -> V  ',8a1)
			    else
			    ilum=3
c			    print 1006,icnt,identifier,otype,obafg,subtype,lum
1006			    format(i6,2x,a26,a4,2x,', Luminosity -> III',8a1)
			    endif
			  endif
c			  
			  if(ibump.ge.1.and.(lum(1+ibump).ne.lums(1).and.lum(1+ibump).ne.lums(2)))then
			    iflag=2
			    if(namesv.le.4)then
			    ilum=5
c			    print 1005,icnt,identifier,otype,obafg,subtype,lum
			    else
			    ilum=3
c			    print 1006,icnt,identifier,otype,obafg,subtype,lum
			    endif
			  endif
c			  
			  if(lum(1+ibump).eq.lums(2))ilum=5
			  if(lum(1+ibump).eq.lums(1).and.lum(2+ibump).eq.lums(2))ilum=4
			  if(lum(1+ibump).eq.lums(1).and.lum(2+ibump).eq.lums(1))ilum=2
			  if(lum(1+ibump).eq.lums(1).and.lum(2+ibump).eq.lums(1).and.lum(3+ibump).eq.lums(1))ilum=3
			  if(lum(1+ibump).eq.lums(1).and.ilum.eq.0)ilum=1			
c
c			print 1007, icnt,identifier,obafg,subtype,lum,names(namesv),asubtype,ilum
1007			format(i6,2x,a26,14x,8a1,3x,a1,f5.1,i3)
97			continue
			if(iflag.eq.1)iflagcnt=iflagcnt+1
			if(iflag.eq.2)ilumcnt=ilumcnt+1
			if(iflag.eq.3)iscnt=iscnt+1
c			if(iflag.ge.1)go to 98			!only include stars without assignments
c			if(iflag.lt.1.or.iflag.eq.2)go to 98	!only write output file for some flagged stars
			jarg=ilum
			iarg=1+10*(namesv-1)+ifix(asubtype)
			if(jarg.lt.1.or.jarg.gt.5)print *,'Bad Luminosity = ',ilum
			if(iarg.lt.1.or.iarg.gt.70)print *,'Bad Spectral type = ',namesv, asubtype
			amagsmei=amagv-alog10(mv2smei(iarg,jarg))/alog10(2.5)
c
c	Print out retained clusters, "C,S,N,&R" stars, and Wolf-Rayet stars.
c			
			if(icluster.eq.1.and.amagsmei.gt.6.0)then
				icluscnt=icluscnt+1
c				print 1008,icnt,identifier,otype,amagsmei
1008				format(i6,2x,a26,a4,2x,',  Cluster discarded, SMEI mag = ',f5.2)
				go to 98
			endif
c			if(icluster.eq.1)print 1009,icnt,identifier,otype,amagsmei
1009				format(i6,2x,a26,a4,2x,',  Cluster  retained, SMEI mag = ',f5.2)
			if(iflag.eq.3)then
				if(amagsmei.gt.6.0)then
				is=is+1
				go to 98
				else
c				        print 1010,icnt,identifier,otype,obafg,subtype,lum,amagsmei
1010				format(i6,2x,a26,a4,2x,', ',8a1,', retained, SMEI mag = ',f5.2)
				endif
			endif				
			if(iflag.eq.4)then
				iwolf=iwolf+1
				if(amagsmei.gt.6.0)then
					iw=iw+1
					go to 98
				else
c					print 1010,icnt,identifier,otype,obafg,subtype,lum,amagsmei
				endif
			endif				
c
			if(amagsmei.gt.6.0)go to 98		!reject anything fainter than Msmei = 6.0
c			
			call double(0,identifier,iout)		!reject multiple entries & selected clusters
			if(iout.eq.1)then
				ioutcnt=ioutcnt+1
				go to 98
			endif
c
c	Entries to be output arrive here...
c	Flag M7 - M9 troublemakers
c
c			if(iarg.ge.68)print 1011,identifier,obafg,subtype,lum,otype,irahour,iramin,rasec,idecdeg,
c     &			idecmin,idecsec,amagb,amagv,amagsmei
1011	format(a26,8a1,2x,a4,2x,i2,1x,i2,1x,f5.2,2x,i3,1x,i2,1x,i2,8x,f5.2,6x,f5.2,6x,f5.2)     
c
c	Flag eclipsing binaries
c
c			if(otype.eq.eclipse(1))print 1012,icnt,identifier,otype,obafg,subtype,lum,amagsmei			
c			if(otype.eq.eclipse(2))print 1012,icnt,identifier,otype,obafg,subtype,lum,amagsmei			
c			if(otype.eq.eclipse(3))print 1012,icnt,identifier,otype,obafg,subtype,lum,amagsmei			
c			if(otype.eq.eclipse(4))print 1012,icnt,identifier,otype,obafg,subtype,lum,amagsmei			
1012				format(i6,2x,a26,a4,2x,', ',8a1,', eclipsing, SMEI mag =',f5.2)
			itype(iarg,jarg)=itype(iarg,jarg)+1
			ikeep=ikeep+1
		write(11,11)obafg,subtype,lum,otype,irahour,iramin,rasec,idecdeg,idecmin,idecsec,amagb,amagv,amagsmei,
     &			iarg,jarg,identifier			
11		 	format(8a1,2x,a4,2x,i2,1x,i2,1x,f5.2,2x,i3,1x,i2,1x,i2,8x,f5.2,6x,f5.2,6x,f5.2,2i4,2x,a26)
98			continue
		enddo
99		continue
c		print 15,(heads(10)),icnt
15		format(' Finished'a50,7x,i6,' stars')
	   enddo
	   irahour=0
	   iramin=0
	   rasec=0.0
	   write(11,11)obafg,subtype,lum,otype,irahour,iramin,rasec,idecdeg,idecmin,idecsec,amagb,amagv,amagsmei
	   close(11)
	   close(12)
	enddo
	print *, icount,' stars considered, ',iflagcnt,' subtypes assigned, '
	print *, ilumcnt,' luminosities assigned,',ikeep,' Retained'
	print *, iscnt,' C S N R-type stars found:',is,' rejected, too faint;'
	print *, iwolf,' Wolf-Rayet stars found:',iw,' rejected, too faint;'
	print *,iibad,' Bad names found'
	print *, icluscnt,' Clusters fainter than 6th mag: discarded'
	print *, ioutcnt, ' Multiple entries or resolved clusters > 6th, deleted'
	print 1020, (j-1,(itype(j,i),i=1,5),j=1,70)
1020	format(i3,10x,5i5)
	open(12,file='sclass.dat')
	write(12,1020) (j-1,(itype(j,i),i=1,5),j=1,70)
c
	stop 'reached end'
	end
***********************************************************************
	subroutine double(ifirst,identifier,iout)
	character*26 identifier,badnames(133)
	integer*4 i,ifirst,iout,imin
	save badnames,imin
c
c	This program removes, by brute force, all identifiers tagged as double/multiple entries
c	and three clusters brighter than 6th that are large enough to be resolved by SMEI
c
	if(ifirst.eq.1)then
	imin=1
	badnames(1)  ='GJ 335A                   '
	badnames(2)  ='HD 223385 A               '
	badnames(3)  ='HD 223385 B               '
	badnames(4)  ='ADS 1507 AB               '
	badnames(5)  ='HD 12447                  '
	badnames(6)  ='HD 12446                  '
	badnames(7)  ='HD 18519 J                '
	badnames(8)  ='HD 34029                  '
	badnames(9)  ='HD 36861 J                '
	badnames(10) ='HD 42126 J                '
	badnames(11) ='HD 60179                  '
	badnames(12) ='HD 60178 J                '
	badnames(13) ='CCDM J08123+1738AB        '
	badnames(14) ='HD 68257                  '
	badnames(15) ='M 44                      '
	badnames(16) ='HD 76943                  '
	badnames(17) ='CCDM J10199+1951AB        '
	badnames(18) ='HD 98231                  '
	badnames(19) ='HD 98230                  '
	badnames(20) ='HD 114378                 '
	badnames(21) ='HD 114379                 '
	badnames(22) ='HD 116656                 '
	badnames(23) ='HD 116657                 '
	badnames(24) ='HD 119425 A               '
	badnames(25) ='HD 124674 J               '
	badnames(26) ='HD 129247                 '
	badnames(27) ='HD 129246                 '
	badnames(28) ='CCDM J14449+2704AB        '
	badnames(29) ='HD 131156 A               '
	badnames(30) ='GJ 575A                   '
	badnames(31) ='HD 137107 J               '
	badnames(32) ='HD 138917 J               '
	badnames(33) ='HD 138917                 '
	badnames(34) ='HD 139891 J               '
	badnames(35) ='HD 146361 J               '
	badnames(36) ='HD 150117 J               '
	badnames(37) ='HD 154906                 '
	badnames(38) ='HD 154905                 '
	badnames(39) ='HD 156014 J               '
	badnames(40) ='HD 157779                 '
	badnames(41) ='HD 157778                 '
	badnames(42) ='GJ 702B                   '
	badnames(43) ='GJ 702A                   '
	badnames(44) ='GJ 704A                   '
	badnames(45) ='IC 4756                   '
	badnames(46) ='HD 173582                 '
	badnames(47) ='HD 173608                 '
	badnames(48) ='HD 173607                 '
	badnames(49) ='GJ 738A                   '
	badnames(50) ='ADS 14636 AB              '
	badnames(51) ='GJ 822A                   '
	badnames(52) ='GJ 822B                   '
	badnames(53) ='CCDM J21442+2845AB        '
	badnames(54) ='HD 213052                 '
	badnames(55) ='HD 213051                 '
	badnames(56) ='HD 3196 A                 '
	badnames(57) ='HD 10360 J                '
	badnames(58) ='HD 18623                  '
	badnames(59) ='HD 18622                  '
	badnames(60) ='HD 24071 J                '
	badnames(61) ='CCDM J03486-3737AB        '
	badnames(62) ='HD 24555                  '
	badnames(63) ='HD 31203                  '
	badnames(64) ='HD 37742 J                '
	badnames(65) ='HD 37743                  '
	badnames(66) ='CCDM J05445-2226AB        '
	badnames(67) ='HD 45725 J                '
	badnames(68) ='CCDM J06288-0702AC        '
	badnames(69) ='HD 45726 J                '
	badnames(70) ='CCDM J06386-4813AB        '
	badnames(71) ='CCDM J07292-4318AB        '
	badnames(72) ='HD 60584                  '
	badnames(73) ='HD 60585                  '
	badnames(74) ='CCDM J07388-2648AB        '
	badnames(75) ='CCDM J08112-1256AE        '
	badnames(76) ='CCDM J08255-5144AB        '
	badnames(77) ='CCDM J08295-4443AB        '
	badnames(78) ='GJ 314A                   '
	badnames(79) ='CCDM J08397-2934AB        '
	badnames(80) ='HD 82434                  '
	badnames(81) ='GJ 366.1B                 '
	badnames(82) ='GJ 366.1A                 '
	badnames(83) ='HD 91355                  '
	badnames(84) ='NGC 3532                  '
	badnames(85) ='CCDM J11323-2916AB        '
	badnames(86) ='HD 100286                 '
	badnames(87) ='HD 100287                 '
	badnames(88) ='HD 110379 J               '
	badnames(89) ='HD 112092                 '
	badnames(90) ='CCDM J13532-3156AB        '
	badnames(91) ='HD 130559                 '
	badnames(92) ='HD 131977                 '
	badnames(93) ='HD 133243                 '
	badnames(94) ='CCDM J15051-4703AB        '
	badnames(95) ='HD 136415                 '
	badnames(96) ='CCDM J15511-5503AB        '
	badnames(97) ='CCDM J16001-3824AB        '
	badnames(98) ='HD 144070                 '
	badnames(99) ='HD 144069                 '
	badnames(100)='CCDM J16054-1948AC        '
	badnames(101)='HD 144218                 '
	badnames(102)='CCDM J16212-2536AB        '
	badnames(103)='HD 147723                 '
	badnames(104)='HD 147933                 '
	badnames(105)='CCDM J16272-4733AB        '
	badnames(106)='HD 148478 J               '
	badnames(107)='HD 155885                 '
	badnames(108)='HD 155886                 '
	badnames(109)='CCDM J17591-3015AB        '
	badnames(110)='HD 164765                 '
	badnames(111)='HD 165190                 '
	badnames(112)='HD 165189                 '
	badnames(113)='HD 177474                 '
	badnames(114)='CCDM J19226-4428AB        '
	badnames(115)='CCDM J19278-5420AB        '
	badnames(116)='CCDM J19299-2659AB        '
	badnames(117)='HD 200496 J               '
	badnames(118)='CCDM J22143-2104AB        '
	badnames(119)='CCDM J22183-5338AB        '
	badnames(120)='HD 213052                 '
	badnames(121)='HD 213051                 '
	badnames(122)='HD 219834                 '
	badnames(123)='* 107 Aqr                 '
	badnames(124)='CCDM J00316-6258AC        '
	badnames(125)='CCDM J03294-6256AB        '
	badnames(126)='HD 85123 J                '
	badnames(127)='HD 99103                  '
	badnames(128)='HD 99103 J                '
	badnames(129)='CCDM J11234-6457AB        '
	badnames(130)='HD 101379        '
	badnames(131)='HD 108248                 '
	badnames(132)='HD 108249                 '
	badnames(133)='HD 128620 J               '
	else
		iout=0
		do i=imin,133
		if(identifier.eq.badnames(i))then
			iout=1
			if(i.ne.imin)print *,i-1,' missed or out of order'
			imin=i+1
c			if(i.ne.15.and.i.ne.44.and.i.ne.81)print *,i,identifier,' discarded, multiple entry'
c			if(i.eq.15.or.i.eq.44.or.i.eq.81)print *,i,identifier,' discarded, resolved cluster'
			return
		endif
		enddo
	endif
c	
	return
	end