pro spectro,fichier
lcar=strlen(fichier)
par=strmid(fichier,1,1)
;type=1 pour B// 2 pour V// 0 pour intensite
type=0
if par eq 'm' then type=1
if par eq 'v' then type=2
mois=strmid(fichier,4,2)
jour=strmid(fichier,6,2)
mois=fix(mois)
jour=fix(jour)
if (mois gt 1) then jour=jour+fix(30.42*float(mois-1))
fichier1=fichier
if type ne 0 then fichier1='td'+strmid(fichier,2,lcar-2)
read_gif,fichier1,tab
dim=size(tab)
im=dim(1)
jm=dim(2)
;
seuil=25
;
jj=jm/4
res=tab(*,jj)
ii=im/2
mar=20
for i=mar,ii do begin
  if res(i) gt seuil then goto, lab1
endfor
lab1: i1=i
for i=im-1-mar,im-ii,-1 do begin
  if res(i) gt seuil then goto, lab2
endfor
lab2: i2=i
j1=jj
j2=j1
ii=im/2
res=tab(ii,*)
jj=jm/2
i3=ii
for j=jm-1-mar,jm-jj,-1 do begin
  if res(j) gt seuil then goto, lab3
endfor
lab3: j3=j
haut1=(i1^2-i2^2+j1^2-j2^2)*(i1-i3)
haut2=(i1^2-i3^2+j1^2-j3^2)*(i1-i2)
bas=(j1-j2)*(i1-i3)-(j1-j3)*(i1-i2)
j0=0.5*(haut1-haut2)/bas
i0=(j1+j2-2*j0)*(j1-j2)/(i1-i2)
i0=0.5*(i1+i2+i0)
r=sqrt((i1-i0)^2+(j1-j0)^2)
;print,i1,i2,i3
;print,j1,j2,j3
;print,i0,j0,r
jj=3*jm/4
res=tab(*,jj)
ii=im/2
for i=mar,ii do begin
  if res(i) gt seuil then goto, lab4
endfor
lab4: i4=i
for i=im-1-mar,im-ii,-1 do begin
  if res(i) gt seuil then goto, lab5
endfor
lab5: i5=i
j4=jj
j5=j4
ii=im/2
res=tab(ii,*)
jj=jm/2
i6=ii
for j=mar,jj do begin
  if res(j) gt seuil then goto, lab6
endfor
lab6: j6=j
haut1=(i4^2-i5^2+j4^2-j5^2)*(i4-i6)
haut2=(i4^2-i6^2+j4^2-j6^2)*(i4-i5)
bas=(j4-j5)*(i4-i6)-(j4-j6)*(i4-i5)
j00=0.5*(haut1-haut2)/bas
i00=(j4+j5-2*j00)*(j4-j5)/(i4-i5)
i00=0.5*(i4+i5+i00)
rr=sqrt((i4-i00)^2+(j4-j00)^2)
;print,i4,i5,i6
;print,j4,j5,j6
;print,i00,j00,rr
i0=(i0+i00)/2.
j0=(j0+j00)/2.
r=(r+rr)/2.
ii0=fix(i0)
jj0=fix(j0)
if type ne 0 then read_gif,fichier,tab
tab(i1,j1)=255
tab(i2,j2)=255
tab(i3,j3)=255
tab(i4,j4)=255
tab(i5,j5)=255
tab(i6,j6)=255
tab(ii0,jj0)=0
tab(ii0-1,jj0)=0
tab(ii0+1,jj0)=0
tab(ii0,jj0-1)=0
tab(ii0,jj0+1)=0
window,0,xsize=im,ysize=jm
tv,tab
if type eq 0 then loadct,0
if type ne 0 then lut,'vit2.lut'
plot,[0,im-1],[0,jm-1],xmargin=[0,0],ymargin=[0,0],/nodata,/device,xstyle=1,ystyle=1,/noerase
xx=fltarr(181)
yy=fltarr(181)
for ind=-90,90 do begin
xx(ind+90)=i0+cos(float(ind)*3.1415926/180.)*r
yy(ind+90)=j0+sin(float(ind)*3.1415926/180.)*r
endfor
oplot,xx,yy
for ind=90,270 do begin
xx(ind-90)=i0+cos(float(ind)*3.1415926/180.)*r
yy(ind-90)=j0+sin(float(ind)*3.1415926/180.)*r
endfor
oplot,xx,yy
;
b0=-7.25*cos(2.*3.1415926*float(jour-65)/365.)
b0=b0*3.1415926/180.
sb=sin(b0)
cb=cos(b0)
p=0.
p=p*3.1415926/180.
sp=sin(p)
cp=cos(p)
;
for ind=1,17 do begin
la=float(ind-9)*10.*3.1415926/180.
sla=sin(la)
cla=cos(la)
for jnd=-90,90 do begin
lg=float(jnd)*3.1415926/180.
slg=sin(lg)
clg=cos(lg)
xx(jnd+90)=i0+r*(cp*cla*slg-sp*sla*cb+sp*cla*clg*sb)
yy(jnd+90)=j0+r*(sp*cla*slg+cp*sla*cb-cp*cla*clg*sb)
endfor
oplot,xx,yy
endfor
for ind=1,17 do begin
lg=float(ind-9)*10.*3.1415926/180.
slg=sin(lg)
clg=cos(lg)
for jnd=-90,90 do begin
la=float(jnd)*3.1415926/180.
sla=sin(la)
cla=cos(la)
xx(jnd+90)=i0+r*(cp*cla*slg-sp*sla*cb+sp*cla*clg*sb)
yy(jnd+90)=j0+r*(sp*cla*slg+cp*sla*cb-cp*cla*clg*sb)
endfor
oplot,xx,yy
endfor
;
on_error,2              ;Return to caller if an error occurs
print,'Press left or center mouse button for new output line."
print,'... right mouse button to exit.'
s = size(tab)
if s(0) ne 2 then message, 'Image parameter not 2d.'
s(1) = s(1)-1		;To n-1
s(2) = s(2)-1
!err=0
if s(s(0)+1) ge 4 then form = 'F' else form = 'I'
cr = string("15b)	;this codes a newline
;
while !err ne 4 do begin
	tvrdc,i,j,2,/dev
	if (!err and 3) ne 0 then begin	;New line?
	   print,form="($,a)",string("12b)
	   while (!err ne 0) do begin wait,.1 & tvrdc,i,j,0,/dev & end
	endif
        ii=i
        jj=j
	i=float(i)
	j=float(j)
        dis=sqrt((i-i0)^2+(j-j0)^2)
        if dis lt r then begin
        xcor=(i-i0)*960.0/r
        ycor=(j-j0)*960.0/r
	phi0=asin((j-j0)/r)
	teta0=asin((i-i0)/(r*cos(phi0)))
	ip=100
	db=b0/float(ip)
	b=0.
	phi=phi0
	teta=teta0
	for ind=1,ip do begin
	haut=sin(b)*sin(phi)+cos(phi)*cos(teta)*cos(b)
	bas=cos(phi)*cos(b)+cos(teta)*sin(b)*sin(phi)+sin(b)*sin(phi)*(sin(teta))^2/cos(teta)
	dphi=db*haut/bas
	dteta=dphi*sin(teta)*sin(phi)/(cos(teta)*cos(phi))
	phi=phi+dphi
	teta=teta+dteta
	b=b+db
	endfor
	phi0=phi0*180./3.1415926
	teta0=teta0*180./3.1415926
	phi=phi*180./3.1415926
	teta=teta*180./3.1415926
;	print,phi0,phi,teta0,teta
        val=tab(ii,jj)
        form='i'
        if type eq 2 then begin
          val=(float(tab(ii,jj))-127.)*0.05
          form='f5.2'
        end
        if type eq 1 then begin
          val1=(float(tab(ii,jj))-127.)*0.03
          if val1 ge 0. then val=10.^val1
          if val1 lt 0. then val=-10.^abs(val1)
          form='f6.0'
        end
	form="($,'latitude=',f5.1,', longitude=',f5.1,', x=',f5.0,', y=',f5.0,', value=',"+form+",a)"
	print,form = form, phi,teta,xcor,ycor,val,cr
        end
endwhile
print,form="(/)"
end