
C-----------------------------------------------------------------------------
C---------------- the Follow is SUBROUTINES in Main Program   ----------------
C-----------------------------------------------------------------------------


C------FSub.3
cc----The program is used to compare the magnetogtram and the Ha picture.-------
cc--Parameters except for the dimension of the spectra are input from specmap---

	Subroutine spec_mag

        parameter (xx1=90.0,xx2=90.,xx3=95.,yy1=60.,yy2=55.,yy3=55.)
        real nxx(3)/xx1,xx2,xx3/,nyy(3)/yy1,yy2,yy3/
        real bx(170,128),by(170,128),tr(6)
        real flz(170,128),rdat(170,128),adat(170,128)
        character fill*12,file*12,filx*12,fily*12,key
        character cht(8)*5
        integeranw8,fl,sl,rl,bl,sp
        CHARACTER MESG*400
        real gy(512,512),tmp(512,512)
        real cm(8),cp(8),fm(8),fp(8)
	real xmo(20),ymo(20),xmg(20),ymg(20),angmo(10),angmg(10)
        real gm(6),trmv(6)
        REAL smap(208,40),trs(6),trm(6),ggy(208,40)
        INTEGER*2 dispec(250,50)
        integer curx1,curx2,cury1,cury2
	real slev(10),flev(10),rlev(10),blev(10),spo(10),lev3(10),lev5(10)
	real sslev,fflev,rrlev,bblev,sspo,resox,resoy,scalx,scaly


c       data gm/90.,140.,170.,200.,240.,256./
c       data cm/1200.,1500.,1800.,2000.,2400.,2800.,3400.,4000./
	data cht/'0.20','0.40','0.80','1.00','1.20','1.40','1.60','2.00'/
c       data fm/1152.0,1071.0,984.0,890.0,677.0,424.0,160.0,40.0/
c       DATA FM/2880.0,2560.0,2240.,1920.,1280.,640.,160.,40.0/
        data tr/0.0,1.0,0.0,0.0,0.0,1.0/
        data trmv/0.0,1.0,0.0,0.0,0.0,1.0/
        DATA trs/0.0,1.0,0.0,0.0,0.0,1.0/
        DATA trm/0.0,1.0,0.0,0.0,0.0,1.0/
        external r_map !$pragma C (r_map)
        external rw_fit !$pragma C (rw_fit)
        external rspecmap !$pragma C (rspecmap)

        pi=3.1415926

        open(2,file='cur_spec.lev',status='old')
        read(2,*)ii0,jj0,resox,resoy
        read(2,*)fl,fflev,(flev(i),i=1,fl)
        read(2,*)sl,sslev,(slev(i),i=1,sl)
        read(2,*)sp,sspo,(spo(i),i=1,sp)
        read(2,*)rl,rrlev,(rlev(i),i=1,rl)
        read(2,*)bl,bblev,(blev(i),i=1,bl)
        read(2,*)icm,(cm(i),i=1,icm)
        read(2,*)igm,(gm(i),i=1,igm)
        read(2,*)ifm,(fm(i),i=1,ifm)
        read(2,*)curx1,curx2,cury1,cury2
        read(2,*)lx0,lxx,ly0,lyy
        read(2,*)nnx0,nnxx,nny0,nnyy
        read(2,*)mtx0,mtxx,mty0,mtyy
 
        do m=1,icm
           cp(m)=-cm(m)
        end do
        do m=1,ifm
           fp(m)=-fm(m)
        end do
        smdx=0.
        smdy=0.
              do i=1,fl
                 flev(i)=fflev*flev(i)
              end do
              do i=1,sl
                 slev(i)=sslev*slev(i)
              end do
              do i=1,sp
                 spo(i)=sspo*spo(i)
              end do
              do i=1,rl
                 rlev(i)=rrlev*rlev(i)
              end do
              do i=1,bl
                 blev(i)=bblev*blev(i)
              end do

		xhc=104.
		yhc=20.
                xc=256.
                yc=256.
		pvxh=1.4*resox
		pvyh=1.4*resoy
                px=5.23/512.
                py=3.626/512.
                pvx=px*1.4
                pvy=py*1.4
		scalx=resox*208.
		scaly=resoy*40.
		scalx=scalx/60.
		scaly=scaly/60.
	        sx0=(5.23-scalx)*1.4/2.+1.0
		sy0=(3.626-scaly)*1.4/2.+0.8
111	write(*,*)'*******************************************************'
	write(*,*)'***<Menu 4 Comparision with 2D-spectr.& Magnetogram>***'
	write(*,*)'*******************************************************'
	write(*,*)'** 0. go to main menu                                **'
	write(*,*)'** 1. rotate an angle for magnetic map               **'
	write(*,*)'** 2. move the spectrum map                          **'
	write(*,*)'** 3. display or print                               **'
	write(*,*)'** 4. Quit                                           **'	
	write(*,*)'*******************************************************'
	     write(*,*)'----select num:(1-4)'
	     read(*,*)anw8
                tr(2)=cosd(theta)
                tr(3)=-(resoy/resox)*sind(theta)
                tr(1)=xhc-xhc*cosd(theta)+yhc*(resoy/resox)*sind(theta)
                tr(5)=(resox/resoy)*sind(theta)
                tr(6)=cosd(theta)
                tr(4)=yhc-xhc*(resox/resoy)*sind(theta)-yhc*cosd(theta)
c                trm(2)=cosd(theta)
c                trm(3)=-(py/px)*sind(theta)
c                trm(1)=xc-xc*cosd(theta)+yc*(py/px)*sind(theta)
c                trm(5)=(px/py)*sind(theta)
c                trm(6)=cosd(theta)
c                trm(4)=yc-xc*(px/py)*sind(theta)-yc*cosd(theta)
	     IF(ANW8.EQ.0) THEN
	        goto 999
	     ELSE IF(ANW8.EQ.1.OR.ANW8.EQ.2) THEN
                write(*,*)'----Input Magnetic-Image File Name---->'
                read(*,*)fill
                call r_map(fill,tmp)
                do i=1,512
                do j=1,512
                 gy(j,i)=tmp(j,513-i)
                end do
                end do
C     read spectrum image file -----------
                call rspecmap(smap)
                do  j=1,40
                do  i=1,208
                   if(smap(i,j).gt.150.)smap(i,j)=0.
        	end do
        	end do
	        do  j=1,40
	        do  i=1,208
	            dispec(i,j)=int(smap(i,j))
	        end do
	        end do
	        write(*,*)'Do you rotate the map before the machine?<y/n>'
		read(*,*)key
		if(key.eq.'y'.or.key.eq.'Y')then
		 write(*,*)'input rotating angle ==>'
		 read(*,*)theta
                 tr(2)=cosd(theta)
                 tr(3)=-(resoy/resox)*sind(theta)
                 tr(1)=xhc-xhc*cosd(theta)+yhc*(resoy/resox)*sind(theta)
                 tr(5)=(resox/resoy)*sind(theta)
                 tr(6)=cosd(theta)
                 tr(4)=yhc-xhc*(resox/resoy)*sind(theta)-yhc*cosd(theta)
		end if
		write(*,*)'Do you move the map before the machine?<y/n>'
		read(*,*)key
		if(key.eq.'y'.or.key.eq.'Y')then
		 write(*,*)'Input dx and dy'
		 read(*,*)smdx,smdy
		 sx0=sx0+smdx
		 sy0=sy0+smdy
		end if
              call fmaxmin(gy,512,512,white,black)
   	      call pgbegin(0,'?',1,1)
              call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
	      call pgbox('bcs',0.0,0,'bcs',0.0,0)
	      call pgwindow(1.,512.,1.,512.)
	      CALL PGGRAY(gy,512,512,1,512,1,512,WHITE,BLACK,TRm)
	      call fmaxmin(smap,208,40,black,white)
	      call pgvsize(sx0,1.4*scalx+sx0,sy0,1.4*scaly+sy0)
	      call pgwindow(1.,208.,1.,40.)
	        write(*,*)'-------Contour map or Gray map?<C/G>'
	        read(*,*)key
	        if(key.eq.'g') then
	         call pggray(dispec,208,40,1,208,1,40,white,black,tr,2)
	        else if(key.eq.'c') then
	          call pgsci(3)
	          call pgslw(3)
	 	  call pgcont(smap,208,40,1,208,1,40,spo,-sp,tr) 
	          call pgslw(1)
	        end if
                call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
                call pgwindow(1.,512.,1.,512.)
C-----------rotate  the magnetic map --------
	      if(anw8.eq.1) then
	        write(*,*)'Do you rotate the magnetic map?<Y/N>'
	        read(*,*)key
	        if(key.eq.'y') then
C-------- select contrast points
                  do k=1,nch
                   xmo(k)=0.0
                   ymo(k)=0.0
                   xmg(k)=0.0
                   ymg(k)=0.0
                end do
                nch=0
	    write(*,*)'---Select position  points on  spec-image map (BLUE)--'
                call pgsci(5)
                call pgolin(20,nch,xmo,ymo,18)
	   write(*,*)'---Select position  points on magnect-image (YELLOW-RED)--'
                nch=0
                call pgsci(4)
                call pgolin(20,nch,xmg,ymg,17)
                write(*,*)'Select point number as contrast:-->',nch
	        ang=0.0
	        do k=1,nch/2
	           angmo(k)=atan2d(ymo(2*k)-ymo(2*k-1),xmo(2*k)-xmo(2*k-1))
	           angmg(k)=atan2d(ymg(2*k)-ymg(2*k-1),xmg(2*k)-xmg(2*k-1)) 
	           ang=ang+(angmo(k)-angmg(k))
	        end do
	        theta=theta+ang/(nch/2)
	        write(*,*)'rotate theta==>',theta,'this time ANG=>',ang*2/nch
	      end if
	      call pgend
	      goto 111 
	    else if(anw8.eq.2) then
 	       	write(*,*)'Do you move the position?<y/n>'
    	        read(*,*)key
	        if(key.eq.'y') then	
C-------- select contrast points
		do k=1,nch
		   xmo(k)=0.0
	 	   ymo(k)=0.0
		   xmg(k)=0.0
		   ymg(k)=0.0
		end do
	        nch=0 
	write(*,*)'---Select position  points on  spec-image map (BLUE)--'
	        call pgsci(5)
                call pgolin(20,nch,xmo,ymo,18)
	write(*,*)'---Select position  points on magnect-image (YELLOW-RED)--'
	        nch=0
	        call pgsci(4)
	        call pgolin(20,nch,xmg,ymg,17)
	        write(*,*)'Select point number as contrast:-->',nch
	        smdx=0.0
	        smdy=0.0
	        do k=1,nch
	           smdx=smdx+(xmg(k)-xmo(k))
		   smdy=smdy+(ymg(k)-ymo(k))
	        end do
	        smdx=pvx*(smdx/nch)
	        smdy=pvy*(smdy/nch)
	        write(*,65)smdx,smdy
65  	        format('move dx= ',f8.3,'  move dy= ',f8.3)
	        sx0=sx0+smdx
	        sy0=sy0+smdy
	        dxall=dxall+smdx
	        dyall=dyall+smdy
	        write(*,*)'sum dx=',dxall,' sum dy=',dyall
	     else
		write(*,*)'BEGIN RITING!!!'
	   write(2,2222)theta
	   write(2,2223)smdx,smdy
2222	   format('last theta ==',f8.3)
2223	   format('last dx ==',f8.3,'last dy ==',f8.3)
		write(*,*)'WRITE OK!!!'
	   close(2)
	    end if
            call pgend	
	    goto 111
	   end if
	  ELSE IF(ANW8.EQ.3) THEN
	     write(*,*)'-----Input Longitude Field File:[*.bzo]--->'
             read(*,*)file
             call rw_fit(file,flz,mesg,'r')
	     call fmaxmin(flz,170,128,black,white)
	     call rspecmap(smap)
c	     call fmaxmin(smap,208,40,specb,specw)
        do  j=1,40
        do  i=1,208
           if(smap(i,j).ge.150.)smap(i,j)=0.
	end do
           smap(208,j)=smap(207,j)
	end do
            xc=170./2.
            yc=128./2.
            px=5.23/170.
            py=3.626/128.
            pvx=5.23*1.4/170.
            pvy=3.626*1.4/128.
	   write(*,*)'----change rotate angle or moved value?<y/n>'
	   read(*,*)key
	   if(key.eq.'y') then
	     write(*,*)'rotate angle==>',theta
	     write(*,*)'moved dx=>',dxall,'moved dy=>',dyall
	     write(*,*)'--input the rotate angle (degree):'
	     read(*,*) theta
	     write(*,*)'--input  the change value  of spectr-map position:(dx,dy)'
	     read(*,*) smdx,smdy
	   end if
            sx0=1.0+smdx+(5.23-scalx)*1.4/2.
            sy0=0.8+smdy+(3.626-scaly)*1.4/2.
                tr(2)=cosd(theta)
                tr(3)=-(resoy/resox)*sind(theta)
                tr(1)=xhc-xhc*cosd(theta)+yhc*(resoy/resox)*sind(theta)
                tr(5)=(resox/resoy)*sind(theta)
                tr(6)=cosd(theta)
                tr(4)=yhc-xhc*(resox/resoy)*sind(theta)-yhc*cosd(theta)
              call pgbegin(0,'?',1,1)
              call pgvsize(sx0,1.4*scalx+sx0,sy0,1.4*scaly+sy0)
              call pgbox(' ',0.0,0,' ',0.0,0)
              call pgwindow(1.,208.,1.,40.)
	    write(*,*)'----Draw Ha-image or Velocity-image <H/V>?'
	    read(*,*)key
	    if(key.eq.'v') then
	      ll3=bl
	      ll5=rl
	      do i=1,bl
	        lev3(i)=blev(i)
	      end do
	      do i=1,rl
	        lev5(i)=rlev(i)
	      end do
	      do i=1,208
	       do j=1,40
		smap(i,j)=-smap(i,j)
	       end do
	      end do
	    else
	      ll3=sl
	      ll5=fl
	       do i=1,sl
	         lev3(i)=slev(i)
	       end do
	       do i=1,fl
	         lev5(i)=flev(i)
	       end do
	    end if 
	    write(*,*)'Contour or Gray?<c/g>'
	    read(*,*)key
	    if(key.eq.'c'.or.key.eq.'C')then
	      call pgsci(3)
	      call pgslw(2)
	      call pgsls(1)
	      call pgcont(smap,208,40,1,208,1,40,lev5,-ll5,tr)
	      call pgsci(5)
	      call gsls(2)
              call pgcont(smap,208,40,1,208,1,40,lev3,-ll3,tr)
	      call pgslw(1)
	    else if(key.eq.'g'.or.key.eq.'G')then

             write(*,*)'Do you want to change contrast?<y/n>'
             read(*,*)key
             IF(KEY.EQ.'y') THEN
        write(*,*)'input Bright=?[0.0] , Black-White Contrast=[1.0]?'
        READ(*,*)gy0,ratio
                 do i=1,40
                 do j=1,208
                    ggy(j,i)=(smap(j,i)+gy0)/255.
                    ggymin=min(ggy(j,i),ggymin)
                    ggymax=max(ggy(j,i),ggymax)
                 end do
                 end do
                 specw=ggymax*ratio
                 specb=ggymin*ratio
              END IF

             call pgsci(1)
	     call pggray(ggy,208,40,1,208,1,40,specb,specw,tr)
             call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
	      call pgbox('bcs',0.0,0,'bcs',0.0,0) 
              call pgwindow(1.,170.,1.,128.)
cc	    write(*,*)'Input magnetogram size -->(x:1,170; y: 1,128)'
cc	    read(*,*)mx0,mxx,my0,myy
	    write(*,*)'Contour or Gray?<c/g>'
	    read(*,*)key
	    if(key.eq.'c'.or.key.eq.'C')then
c             call pgsci(3)
	      call pgslw(1)
	      call pgsls(1)
	      CALL PGSCI(0)
	      call pgcont(flz,170,128,lx0,lxx,ly0,lyy,fm,-ifm,trm)
c            call pgsci(5)
	      call gsls(2)
	      call pgcont(flz,170,128,lx0,lxx,ly0,lyy,fp,-ifm,trm)
	      call pgslw(1)
	    else if(key.eq.'g'.or.key.eq.'G')then
              call pggray(flz,170,128,lx0,lxx,ly0,lyy,white,black,trm)
	    end if
	    end if
             call pgsci(2)
	      call pgsls(1)
             CALL PGSLW(5)
             CALL PGCONT(flz,170,128,nnx0,nnxx,nny0,nnyy,0.,-1,TRm)

	     WRITE(*,*)'--- Display the Transverse Field?<Y/N>'
	     READ(*,*)KEY
	     IF(KEY.EQ.'y') THEN
	        WRITE(*,*)'--- INPUT [*.bxo]-->'
	        READ(*,*)FILX
	        CALL RW_FIT(FILX,BX,mesg,'r')
	        WRITE(*,*)'--- INPUT [*.byo]-->'
	        READ(*,*)FILY
	        CALL RW_FIT(FILY,BY,mesg,'r')
	        DO I=1,128
	        DO J=1,170
		   RDAT(J,I)=SQRT(BX(J,I)**2+BY(J,I)**2)
		   ADAT(J,I)=ATAN2D(BY(J,I),BX(J,I))
		END DO
		END DO
		CALL FMAXMIN(RDAT,170,128,FMAX,FMIN)
		dqu=1.5*FMAX
	        CALL PGSLW(3)
		CALL MVARROW(RDAT,ADAT,mtx0,mtxx,mty0,mtyy,dqu,trm(1),trm(4))
	    END IF
	     call pgend
	      goto 111
	   ELSE IF(ANW8.EQ.4) THEN
	     stop 
	   END IF
999       close(2) 
	END
	

C---------FSub.4. the Calibration for Stokes Parameters V, Q and U
	subroutine calib(gdat,rdat,adat,dqu,mesg)
	REAL gdat(170,128),rdat(170,128),adat(170,128)
	CHARACTER mesg*400
        real IL(512,512),IU(512,512),IQ(512,512)
        real qdat(170,128),udat(170,128)
	real interp !interpolate function
        character FLL*12,FLQ*12,fLU*12
	character DATE*11,TIMEL*8,SS,key,dew,dsn,temp*9
        character*6 L,AL,BT,BBC,L,XX,YY,FR
	real POSIT(7)
        CHARACTER CDATE*30,SEEING*15,WAVE*25,CFILE*51,CONT*35,SUN*25
        CHARACTER CARR*40,FRAME*30,SUNDISK*40,ORTHO*40
	external r_data !$pragma C (r_data)
	external r_head !$pragma C (r_head)

C-------- Read Data and Do (3 X 4) Smooth: 512 X 512 ===> 170 X 128  ----------
C-------- the Observed Area Range in (5'.23X3'.626) of Solor Disk  ----------
	write(*,*)'--- Input  L File Name ---'
	read(*,*)FLL
	call r_data(FLL,IL)
	call r_head(FLL,date,timel,ss,al,bt,bbc,l,xx,yy,fr,posit)
        call smooth(IL,gdat)
	write(*,*)'Read Q-data and U-data?<Y/N>'
	read(*,*)key
	if(key.eq.'y'.or.key.eq.'Y') then
	  write(*,*)'--- Input  Q File Name ---'
	  read(*,*)FLQ
	  call r_data(FLQ,IQ)
          call smooth(IQ,qdat)
	  write(*,*)'--- Input  U File Name ---'
	  read(*,*)FLU
	  call r_data(FLU,IU)
	  call smooth(IU,udat)
	end if
c	write(*,*)'************** Ok1! **************'

C------  some Parameters about the Recorded Time,the Location and Disk Coord.
C----- of the AR. the East or North below imply Observe Coord.  ('N' mean the
C-----  north extrame of Earth's Rotion Axis.        -------------------------
        sun='Sun Disk Coordinates  : '
        if(posit(2).lt.0.) then
           dew='E'
           n=2
        end if
        if(posit(2).gt.0.) then
           dew='W'
           n=1
        end if
        if(posit(1).lt.0.) then
           m=2
           dsn='S'
        end if
        if(posit(1).gt.0.) then
                m=1
                dsn='N'
        end if
	SUNDISK=sun//dsn//AL(m:5)//' '//dew//BT(n:5)
	ORTHO='Orthogonal Coordinates : '//XX(1:6)//'  '//YY(1:6)
	CARR='Carrington Coordinates : '//BBC(1:5)//'  '//L(1:5)
        n=3
        if(posit(7).ge.1000.0) n=4
        FRAME='Frame Number: '//FR(1:n)
        WAVE='Wave Length : \\gl5324 \\A'
        SEEING='Seeing :  '//SS//' '
        CDATE='Time : '//date//' '//TIMEL
        FLL(1:1)='l'
        write(*,*)'  filename====>',FLL
	temp=FLL(1:9)
	write(*,*)'DATE===> ',cdate
        CFILE='File Name : '//FLL//' '//FLQ//' '//FLU
        CONT='Contents  : Vector Magnetic Field'
C------ put message of data file into MESG*80 -----
	MESG(1:30)='HUAIROU BEIJING ASTRO. OBS. '
	MESG(30:60)=cdate	
	MESG(60:110)=cfile
	MESG(110:150)=sundisk
	MESG(150:190)=ortho
	MESG(190:230)=carr
	MESG(230:260)=frame
	MESG(260:280)=seeing
	MESG(280:320)=cont
	MESG(320:350)=wave
	WRITE(*,*) MESG

	
C--------------    the Selected Frame             ----------------------------
	ny1=1
	ny2=128
	nx1=1
	nx2=170
        X1=FLOAT(NX1)
        X2=FLOAT(NX2)
        Y1=FLOAT(NY1)
        Y2=FLOAT(NY2)

C------  Calculate the Size and Direction of Transverse Field ----------------
	btmax=-1.e6
        do i=ny1,ny2
        do j=nx1,nx2
	     qq=qdat(j,i)
             uu=udat(j,i)
	     rdat(j,i)=sqrt(qq**2+uu**2)
	     if(btmax.lt.rdat(j,i)) btmax=rdat(j,i)
             if(qq.eq.0.0) then
              if(uu.gt.0.0) adat(j,i)=90./2.
              if(uu.lt.0.0) adat(j,i)=3*90./2.
             else 
	       adat(j,i)=90-0.5*atan2d(uu,qq)
             endif
	end do 
	end do
	write(*,*)'no cail-->btmax:',btmax

C_____ two kind of methods about cailbration  about transversional feild____
C------ 1.the old cailbration: Bt=60000*(Q**2+U**2)**(1/2)  (Guass)
C-----  2.the new cailbration: Bt= 9370*(Q**2+U**2)**(1/4)  (Guass)
C--Notice! qq=20000*Q  and uu=20000*U . So cail_old=3.0 and cail_new=93.7/1.414
	write(*,*)'Old (1/2) or New (1/4)? <O/N>'
	read(*,*)key
	if(key.eq.'o') then
	  cail=3.0
	  do i=1,128
	  do j=1,170
	     rdat(j,i)=cail*rdat(j,i)
	  end do
	  end do
	  rmax=cail*btmax
	  write(*,10)rmax
10        format('btmax==>',f8.2,' Guass')
	else if(key.eq.'n') then
	  cail=93.7/1.414
	  do i=1,128
	  do j=1,170
	     rdat(j,i)=cail*sqrt(rdat(j,i))
	  end do
	  end do
	  rmax=cail*sqrt(btmax)
	  write(*,10)rmax
	end if

	dqu=1.2*rmax

C------------------cailbration about the longitudinal field ---------
C-----  Using the Curve of cailbration by Interpolation method -----
	gmax=-1.e6
	gmin=1.e6
c	WRITE(*,*)'---INPUT CAILBRATION OF LONGITUDINAL FIELD==>[2000]:'
c	READ(*,*)CAIL
c	CAIL=1000./cail	
	cail=0.5
	do i=1,128
	do j=1,170
	   gdat(j,i)=cail*gdat(j,i)
           if(gdat(j,i).lt.0.0) then
              gdat(j,i)=-interp(-gdat(j,i))
           else
              gdat(j,i)=interp(gdat(j,i))
           end if
	   if(gmax.lt.gdat(j,i)) gmax=gdat(j,i)
	   if(gmin.gt.gdat(j,i)) gmin=gdat(j,i)
	end do 
	end do
	write(*,20)gmax,gmin
20      format('bzmax==>',f8.2,' (Guass) bzmin==>',f8.2,' (Guass)')
	END

C------FSub.5. Judge Direction of Transverse Field***POTENTIAL  FIELD METHOD **
	SUBROUTINE poten(bzo,rdat,adat,prdat,padat,tmax)
	REAl bzo(170,128),adat(170,128),rdat(170,128)
	real prdat(170,128),padat(170,128)
	real bzi(128,128),bzr(128,128),bxp(128,128),byp(128,128)
	real bxi(128,128),byi(128,128)
 	integer mx,my
	character key
	pi=3.1415926
	mx=128
	my=128
	do i=1,128
        do j=1,128
	   lj=j+20
         bzi(j,i)=0.0
         bzr(j,i)=bzo(lj,i)
        enddo
        enddo
c
        eta=1.
        call d2ft(7,128,eta,bzr,bzi)
c
        do j=1,my
         do i=1,mx
          ai=float(i)
          aj=float(j)
	  rk=sqrt(sin(2*pi*ai/mx)**2/1.**2+sin(2*pi*aj/my)**2/1.**2)
          bxp(i,j)=-bzi(i,j)*sin(2*pi*ai/mx)/(rk*1.)
          bxi(i,j)=-bzr(i,j)*sin(2*pi*ai/mx)/(rk*1.)
          byp(i,j)=-bzi(i,j)*sin(2*pi*aj/my)/(rk*1.)
          byi(i,j)=-bzr(i,j)*sin(2*pi*aj/my)/(rk*1.)
         enddo
        enddo
c
        eta=-1.0
        call d2ft(7,128,eta,bxp,bxi)
        call d2ft(7,128,eta,byp,byi)
c
        do i=1,my
         ii=my+1-i
        do j=1,mx
         jj=mx+1-j
          bxi(j,i)=bxp(jj,ii)
          byi(j,i)=byp(jj,ii)
        enddo
        enddo
c
         call btba(bxi,byi,prdat,padat,tmax)

c	ADJUST THE DIRECTION OF THE OBJECTED TRANSVERTION FIELDS
	WRITE(*,*)'Adjust Direction of Transverse Field?<Y/N>'
	READ(*,*)KEY
	IF(KEY.EQ.'y') THEN
	  tnoise=40.
	  do i=1,128
	  do j=1,128
	     lj=j+20
c	     if(rdat(lj,i).gt.tnoise) then
	        ddd=bxi(j,i)*cosd(adat(lj,i))+byi(j,i)*sind(adat(lj,i))
	        if(ddd.lt.0.) adat(lj,i)=adat(lj,i)+180
c	     end if
	  end do
	  end do
	END IF
	return
	END

C------- FSub.6. Judge the Direction with ********THE CHARGE METHOD**********
	SUBROUTINE charge(gdat,rdat,adat,prdat,padat,tmax)
	REAL rdat(170,128),adat(170,128),gdat(170,128)
	real prdat(170,128),padat(170,128)
	lbox=10
	write(*,*)'give the box size for calculation:(1-20)===>[10] '
        read(*,*)lbox
	tmax=0.0
        do i=1+lbox,128-lbox
        do j=1+lbox,170-lbox
          bbx=0.0
          bby=0.0
          do k=i-lbox,i+lbox
           do l=j-lbox,j+lbox
             if(k.eq.i.and.l.eq.j) goto 123
             dx=float(j-l)*725.*1.
             dy=float(i-k)*725.*1.
             bbx=bbx+gdat(l,k)/(dx**2+dy**2)**1.5*dx
             bby=bby+gdat(l,k)/(dx**2+dy**2)**1.5*dy
123         continue
            enddo
           enddo
	  tt=1.0e5*sqrt(bbx*bbx+bby*bby)
	   prdat(j,i)=tt
	   if(tt.gt.tmax) tmax=tt
	  if(bbx.eq.0.) then
	     if(bby.gt.0.) padat(j,i)=90.
	     if(bby.lt.0.) padat(j,i)=270.
	  else
	    padat(j,i)=atan2d(bby,bbx)
	  end if 
          ddd=bbx*cosd(adat(j,i))+bby*sind(adat(j,i))
          if(ddd.lt.0.0) adat(j,i)=180+adat(j,i)
        enddo
        enddo
	return
	end

C--------FSub.7. Judge the Direction with**** THE KRALL MEHTOD ********
        SUBROUTINE krall(gdat,rdat,adat,prdat,padat,tmax)
        REAL gdat(170,128),rdat(170,128),adat(170,128)
	real prdat(170,128),padat(170,128)
        tnoise=40.
	tmax=0.0
        do i=2,127
        do j=2,169
              bzx=(gdat(j+1,i)-gdat(j-1,i))/2./725.
              bzy=(gdat(j,i+1)-gdat(j,i-1))/2./725.
	      bbx=-bzx
	      bby=-bzy
	      tt=2000.*sqrt(bbx*bbx+bby*bby)
	      prdat(j,i)=tt
	     if(tt.gt.tmax) tmax=tt
	     if(bbx.eq.0.) then
	       if(bby.gt.0.) padat(j,i)=90.
	       if(bby.lt.0.) padat(j,i)=270.
	      else
	       padat(j,i)=atan2d(bby,bbx)
	    end if 
             ddd=cosd(adat(j,i))*bzx+sind(adat(j,i))*bzy
             if(ddd.gt.0.0) adat(j,i)=adat(j,i)+180.
        end do
        end do
        return
        END

C------- FSub.8. Judge the Direction with *******THE DIVERGE METHOD **********
        SUBROUTINE diver(gdat,rdat,adat)
        REAL gdat(170,128),rdat(170,128),adat(170,128)
        tnoise=40.
        do i=2,127
        do j=2,169
          if(rdat(j,i).gt.tnoise) then
            bx1=rdat(j+1,i)*cosd(adat(j+1,i))
            bx2=rdat(j-1,i)*cosd(adat(j-1,i))
            bxx=(bx1-bx2)/2.
            by1=rdat(j,i+1)*sind(adat(j,i+1))
            by2=rdat(j,i-1)*sind(adat(j,i-1))
            byy=(by1-by2)/2.
            ddd=gdat(j,i)*(bxx+byy)
            if(ddd.lt.0.) adat(j,i)=adat(j,i)+180.
          end if
        end do
        end do
        return
        END


C-------FSub.9.	---Adjust the Direction for an Irregular Area by Cursor
	SUBROUTINE adcur(rdat,adat,jx1,jx2,iy1,kl)
        real rdat(170,128),adat(170,128)
        integer jx1(kl),jx2(kl),iy1(kl)
c
        Real xpt(50),ypt(50)
        REAL xsc(50),ysc(50),xsk(50),ysk(50)

	angd=20.
        cang=cosd(angd)
c
        write(*,*)'--- Select N-pole points!'
        npt=0
         call pgsci(5)
         call pgolin(50,npt,xpt,ypt,18)
        NSC=NPT
        do n=1,NSC
         xsc(n)=xpt(n)
         ysc(n)=ypt(n)
        enddo
c
c       now epecify the sink !!!

        write(*,*)'--- select S-pole points!'
        npt=0
         call pgsci(4)
         call pgolin(50,npt,xpt,ypt,17)
        NSK=NPT
        do n=1,NSK
         xsk(n)=xpt(n)
         ysk(n)=ypt(n)
        enddo
C
c for new window to determine the direction
c
c       write(*,*)'************** Ok6! **************'
 
        do m=1,kl
          i=iy1(m)
          jj1=jx1(m)
          jj2=jx2(m)
         do j=jj1,jj2
          IF(RDAT(J,I).GE.TNOISE) THEN
           dscmin=10000000.
           DO N=1,nsc
            dcx=float(j)-xsc(n)
            dcy=float(i)-ysc(n)
            ddc=sqrt(dcx**2+dcy**2)
            if(ddc.lt.dscmin) then
             dscmin=ddc
             nscmin=n
            endif
           enddo
c
          dskmin=10000000.
           DO N=1,nsk
            dkx=xsk(n)-float(j)
            dky=ysk(n)-float(i)
            ddk=sqrt(dkx**2+dky**2)
            if(ddk.lt.dskmin) then
             dskmin=ddk
             nskmin=n
            endif
           enddo
          scx=float(j)-xsc(nscmin)
          scy=float(i)-ysc(nscmin)
          scd=sqrt(scx**2+scy**2)
          scc=scx/scd*cosd(adat(j,i))+scy/scd*sind(adat(j,i))
c
          skx=xsk(nskmin)-float(j)
          sky=ysk(nskmin)-float(i)
          skd=sqrt(skx**2+sky**2)
          skc=skx/skd*cosd(adat(j,i))+sky/skd*sind(adat(j,i))
c
          if(scc.lt.0.0.and.skc.lt.0.0) then
 	      adat(j,i)=180.+adat(j,i)
	  end if
           if(scc.lt.0.0.and.skc.gt.0.0) then
            if(abs(scc).gt.cang) adat(j,i)=180.+adat(j,i)
           endif
           if(skc.lt.0.0.and.scc.gt.0.0) then
            if(abs(skc).gt.cang) adat(j,i)=180.+adat(j,i)
           endif
          endif
         enddo
        enddo
c       write(*,*)'************** Ok7! **************'
        RETURN
	END

C------FSub.10. Adjust the Direction with Cursor for Box Area
	SUBROUTINE box_adcur(rdat,adat,cbx,cby,baa,lm)
	real xpt(2),ypt(2),cbx(50),cby(50),baa(50)
	real rdat(170,128),adat(170,128)
	character key
	lm=0
	do k=1,50
	  write(*,*)' begin.......'
	  npt=0
	  call pgsci(1)
	  call pglcur(2,npt,xpt,ypt)
	  write(*,*)'xpt,ypt===>',xpt,ypt
	  dx=xpt(2)-xpt(1)
	  dy=ypt(2)-ypt(1)
	  if(dx.eq.0.0) then
	    if(dy.gt.0.0) baa(k)=90.
	    if(dy.lt.0.0) baa(k)=270.
	  end if
	  baa(k)=atan2d(dy,dx)
	  cbx(k)=10*int(xpt(1)/10)+5
	  cby(k)=10*int(ypt(1)/10)+5
	write(*,*)'k==> cbx==> cby==> baa==>',k,cbx(k),cby(k),baa(k)
	  lm=k
	  write(*,*)'DO YOU WANT TO SELECT NEXT SMALL BOX?<Y/N>'
	  read(*,*)key
	  if(key.eq.'n') go to 111
	end do
111	do n=1,lm
	   lx0=int(cbx(n))-5
	   ly0=int(cby(n))-5
	    do  i=0,10
	    do  j=0,10
		li=ly0+i
		lj=lx0+j
		dxx=cosd(adat(lj,li)-baa(n))
		if(dxx.lt.0.0) then
	           adat(lj,li)=adat(lj,li)+180.
c		   write(*,*)'lj==> li==> ',lj,li
		end if
		if(adat(lj,li).gt.360.) adat(lj,li)=adat(lj,li)-360.
	    end do
	    end do
	end do
	return
	END


C------- FSub.11. Modified Method: djust Direction with continouse Magnetic line
	SUBROUTINE new_magline(rdat,adat,padat,bb,nbx,nby,ll)
        real rdat(170,128),adat(170,128),padat(170,128)
        real nbx(150),nby(150),bb(150)

          amax=-1.e6
          amin=1.e6
          do i=1,120
          do j=1,170
             if(amax.lt.rdat(j,i)) amax=rdat(j,i)
             if(amin.gt.rdat(j,i)) amin=rdat(j,i)
            if(adat(j,i).gt.360.) adat(j,i)=adat(j,i)-360.
            if(adat(j,i).lt.0.)   adat(j,i)=adat(j,i)+360.
            if(padat(j,i).lt.0.) padat(j,i)=padat(j,i)+360.
            if(padat(j,i).gt.360.) padat(j,i)=padat(j,i)-360.
          end do   
          end do
        write(*,*)'Transverse field strength: Min-->Max'
	write(*,*)amin,amax
	write(*,*)'the Points counted with the least value===>[>200gauss]'
	read(*,*)tnoise
	write(*,*)'*** the Selected Area with Min SHEAR ANGLE:(0-90)==>[30deg]'
        read(*,*)ang
	write(*,*)'*** the least Percent of points counted:(N%)==>[50]'
        read(*,*)nm
	ll=0
        do i=1,120,10
        do j=1,128,10
           daa=0.0 
           sumad=0.0
           kk=0    
	   sumxx=0.0
	   sumyy=0.0
	   sumra=0.0
           kki=(i+9)/10
           kkj=(j+9)/10
           do m=1,10
           do n=1,10
             ii=i-1+m
             jj=j-1+n
             if(rdat(jj,ii).gt.tnoise) then
                 aa=cosd(adat(jj,ii)-padat(jj,ii))
c	         if(aa.lt.0.) write(*,*)'no adjust---->(j,i):',kkj,kki
	         daa=daa+acosd(abs(aa))
		 sumxx=sumxx+cosd(adat(jj,ii))
		 sumyy=sumyy+sind(adat(jj,ii))
		 sumra=sumra+rdat(jj,ii)
		 kk=kk+1
             end if
           end do
           end do
           if(kk.ne.0) then
              avda=daa/float(kk)
	      avra=sumra/float(kk)
	      avrr=sqrt(sumxx*sumxx+sumyy*sumyy)/float(kk)
	      avaa=atan2d(sumyy,sumxx)
	      ratio=avrr/avra
            else
              avda=0.
              avaa=0.
            end if
	    if(kk.gt.nm) then
              write(*,*)'Location(j,i):',kkj,kki
	      write(*,*)'Ratio= Sqrt(sum(bx)**2+sum(by)**2)/Sum(Bt(j,i))'
	      write(*,*)ratio
	      write(*,*)'locate(j,i)=',kkj,kki
	      write(*,*)'Avg. oriention:--- Avg. shear angle: ---number counted: '
	      write(*,*)avaa,avda,kk
	    end if
            if(avda.gt.ang.and.kk.gt.nm) then
             ll=ll+1
              nbx(ll)=j+4
              nby(ll)=i+4
              bb(ll)=avaa
             do m=1,10
             do n=1,10
               ii=i-1+m
               jj=j-1+n
               dd=cosd(adat(jj,ii)-avaa)
               if(dd.lt.0.0) adat(jj,ii)=adat(jj,ii)+180.
               if(adat(jj,ii).gt.360.) adat(jj,ii)=adat(jj,ii)-360.
             end do
             end do
           end if
        end do
        end do
        return
        END


C------------------------FSub.12   Flux for Longitude Magnetic ---------------
	Subroutine L_flux(gdat,jx1,jx2,iy,num)
	REAL GDAT(170,128)
	INTEGER JX1(NUM),JX2(NUM),IY(NUM),NUM,COLOR
	character key

100	write(*,*)'------ You select (+) or (-) parity area?<P/N>'
	read(*,*)key	
	if(key.eq.'p') then
	  color=7
	else
	  color=2
	end if

	dx=5.23*60./170.
        dy=3.626*60./128.
	curt=0.0
 
	write(*,*)'----- Give the Min Limit of Integration[40 Guass] ==>'
        read(*,*)vmin
        do m=1,num
          i=iy(m)
          jj1=jx1(m)
          jj2=jx2(m)
          do j=jj1,jj2
           if(color.eq.7) then
             if(gdat(j,i).gt.vmin) then
               curt=curt+gdat(j,i)
               yi=float(i)
               xj=float(j)
               call pgsci(color)
               call pgpoint(1,xj,yi,17)
             end if
           else if(color.eq.2) then
             if(gdat(j,i).lt.-vmin) then
               curt=curt+gdat(j,i)
               yi=float(i)
               xj=float(j)
               call pgsci(color)
               call pgpoint(1,xj,yi,17)
             end if
           end if
          end do
        end do
	flux=curt*dx*dy*(7.26e7)**2
        write(*,5)flux
5	format('magnetic flux---> ',e11.4,' (Guass*cm*cm)')
	return
	end
	

C---------------- FSub.13. Flux of transverse field for a segnment	
	SUBROUTINE T_FLUX(RDAT,ADAT)
	REAL RDAT(170,128),ADAT(170,128)
	real xpt(2),ypt(2)
	character key
	
230	write(*,*)'------ Please select a segnment :left point and right point'
	npt=0
	tflu=0.0
	dx=5.23*60./170.
	k=0
	call pgsci(1)
	call pgolin(2,npt,xpt,ypt,17)
	slnt=-(xpt(2)-xpt(1))/(ypt(2)-ypt(1))
	ang=atand(slnt)
        write(*,*)'verticle line angle==',ang 
	jx1=nint(xpt(1))
	jx2=nint(xpt(2))
	nstep=1
	if(jx1.gt.jx2) then
	   nstep=-1
	end if
	do jx=jx1,jx2,nstep
	   iy=nint(ypt(1)+(jx-jx1)*(-1./slnt))
	   call pgpoint(1,float(jx),float(iy),17)
	   tflu=tflu+rdat(jx,iy)*cosd(adat(jx,iy)-ang)
	   aa=acosd(cosd(adat(jx,iy)-ang))
	   k=k+1
	   write(*,*)'num=   angle==',k,aa
	end do
	sl=dx*7.24e7/abs(cosd(90.-ang))
	xlength=sl*abs(jx2-jx1+1)
	write(*,*)'Segnment length (cm) ===',xlength
	sum=tflu*sl
        write(*,15)sum
15	format('transverse magnetic flux---> ',e11.4,' (Guass*cm)')
        write(*,*)'---- You want to select next segnment?<Y/N>'
        read(*,*)key
        if(key.eq.'y') goto 230
	return
	end


C------------     FSub.14.   Shear angle for transverse field ----------------
	SUBROUTINE SHEAR(adat,rdat,padat,jx1,jx2,iy,num)
	REAL ADAT(170,128),RDAT(170,128),PADAT(170,128)
	INTEGER JX1(NUM),JX2(NUM),IY(NUM),NUM
        real shang(100*100)

	sum=0.0
	swet=0.0
	npos=0
	wet=0.0
	do m=1,num
          i=iy(m)+1
          jj1=jx1(m)
          jj2=jx2(m)
          do j=jj1,jj2
	     val=cosd(adat(j,i)-padat(j,i))
	     ang=acosd(val)
	     sum=sum+ang
	     swet=swet+ang*rdat(j,i)
	     npos=npos+1
	     wet=wet+rdat(j,i)
	     shang(npos)=ang
	     write(*,*)'ang==', npos,ang
	  end do
	end do

	avg=sum/float(npos)
        write(*,*)'avg== ',avg
	squ=0.0
	do k=1,npos
	   squ=squ+(shang(k)-avg)**2
	end do
	delta=sqrt(squ/float(npos))
	write(*,*)'delta== ',delta
	write(*,*)'*************************************************'

        avwet=swet/wet 
        write(*,*)'weight avg== ',avwet 
        squ=0.0 
        do k=1,npos 
           squ=squ+(shang(k)-avwet)**2
        end do
        delta=sqrt(squ/float(npos))
        write(*,*)'delta== ',delta

	return
	end 
	
	
C------FSub.15  annote the label on the contour map of vector field
	Subroutine labcn(gdt,rdt,adt,dqu,mesg)
	REAL gdt(170,128),rdt(170,128),adt(170,128)
	CHARACTER mesg*400

	real gdat(170,128),rdat(170,128),adat(170,128)
        character boxtext(50)*100,key1,key2,key3,key
        integer textlen(2),boxlen(50)
        real xbox(50),ybox(50),simlen(50)
	REAL TR(6),CM(7),CP(7)
	REAL box(4,20)
        integer nbox
 
        DATA TR/0.0,1.0,0.0,0.0,0.0,1.0/
        data cm/-2880.0,-2240.,-1920.,-1280.,-640.,-160.,-40./
	data cp/2880.0,2240.,1920.,1280.,640.,160.,40./
        external rw_fit !$pragma C (rw_fit)
	external get_str !$pragma C (get_str)


        IDIM=170
	JDIM=128
 	I1=1
	I2=170
	J1=1
	J2=128

	DO I=I1,I2
	DO J=J1,J2
	   GDAT(I,J)=GDT(I,J)
	   RDAT(I,J)=RDT(I,J)
	   ADAT(I,J)=ADT(I,J)
	ENDDO
	ENDDO
C--------------    the Selected Frame             ----------------------------
              do i=1,50
                 boxtext(i)=' '
              end do  
             X1=FLOAT(I1)
             X2=FLOAT(I2)
             Y1=FLOAT(J1)
             Y2=FLOAT(J2)

             nbox=0 
  100         write(*,*)'DO YOU WANT TO MARK LEBALS ON MAP?<Y/N>' 
              call get_str(key1,textlen) 
              if(key1.eq.'y'.or.key1.eq.'Y') then
                 write(*,*)'You want to Select a Clip Area?<Y/N>'
                 read(*,*)key
                 if(key.eq.'y'.or.key.eq.'Y') then
		    write(*,*)'Please input (I1,I2)[0-170] (J1,J2)[0-128]:'
                    read(*,*)I1,I2,J1,J2
                 endif

                CALL PGBEGIN(0,'?',1,1)
                call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
                call pgwindow(x1,x2,y1,y2)
                call pgbox('bcs',0.0,0,'bcs',0.0,0)
     	        call narrow(rdat,adat,I1,I2,J1,J2,dqu,box,nbox)
                CALL GSLS(2)
                call pgsci(2)
                CALL PGSLW(1)
                CALL GCONT(gdat,IDIM,JDIM,I1,I2,J1,J2,CM,-7,TR,box,nbox)
                call pgsls(1)
                call pgsci(7)
                CALL PGSLW(1)
                CALL GCONT(gdat,IDIM,JDIM,I1,I2,J1,J2,CP,-7,TR,box,nbox)
                call pgscf(2)
                call pgsch(0.8)
                call pgslw(2)
                call pgsci(1)

                do k=1,nbox
                   write(*,*)boxtext(k),"length==>",textlen(1)
                   call pgtext(box(1,k)+2.,box(3,k)+2.,boxtext(k))
                end do


                 write(*,*)'How many strings do you expect to lebal?'
                 read(*,*)nbox
  200            write(*,*)'label the  string?<y/N>'               
                 call get_str(key2,textlen)
                 if(key2.eq.'y') then  
                    write(*,30)nbox
  30                format('Select string No.[1--',I2,']:')
                    read(*,*)ns
                    xbox(ns)=0.0
                    ybox(ns)=0.0
                    nch=0
		    write(*,*)'Please locate the position you want to mark!'
                    call pgsci(5)
                    call pgolin(1,nch,xbox(ns),ybox(ns),19)
                    call pgtext(xbox(ns)-1.,ybox(ns)-1.,char(48+ns))
                    write(*,*)'Please input the TEXT in box',ns
                    call get_str(boxtext(ns),textlen)
                    boxlen(ns)=textlen(1)
                    n=boxlen(ns)+1
                    boxtext(ns)(n:n)=' '
                    write(*,*)'how many simbols do you input?'
                    read(*,*)simlen(ns)
                    goto 200
                  end if
                  goto 100
              end if
             do k=1,nbox
                box(1,k)=xbox(k)-4.
                box(2,k)=box(1,k)+2.5*simlen(k)
                box(3,k)=ybox(k)-4.
                box(4,k)=box(3,k)+6.
	        write(*,*)box(1,k),box(2,k),box(3,k),box(4,k)
	        write(*,*)boxtext(k)
             end do
	    call pgend

	        CALL PGBEGIN(0,'?',1,1)
                call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
                call pgwindow(x1,x2,y1,y2)
                call pgbox('bcs',0.0,0,'bcs',0.0,0)
                call narrow(rdat,adat,I1,I2,J1,J2,dqu,box,nbox)
                CALL GSLS(2)
                call pgsci(2)
                CALL PGSLW(1)
                CALL GCONT(gdat,IDIM,JDIM,I1,I2,J1,J2,CM,-7,TR,box,nbox)
                call pgsls(1)
                call pgsci(7)
                CALL PGSLW(1)
                CALL GCONT(gdat,IDIM,JDIM,I1,I2,J1,J2,CP,-7,TR,box,nbox)

              call pgscf(2)
              call pgsch(0.8)
              call pgslw(2)
	      call pgsci(1)
              do k=1,nbox
                  write(*,*)boxtext(k),"length==>",textlen(1)
                  call pgtext(box(1,k)+1.5,box(3,k)+2.,boxtext(k))
              end do
              call pgend
                
              write(*,*)'quit or redisplay?<q/r>'
              call get_str(key3,textlen)
              if(key3.eq.'r') then
                goto 100
              end if

	  END

C------FSub.16  anote the labels on mono-image
	Subroutine labgy
	parameter (IDIM=512,JDIM=512)
        real tr(6),box(4,50),simlen(50)
        character fill*12,key1,key2,key3,key
	character boxtext(50)*100,meg*100
        integer nbox,textlen(2),boxlen(50)
	real xbox(50),ybox(50)
        real gy(512,512),tmp(512,512)
	data tr/0.0,1.0,0.0,0.0,0.0,1.0/

        external r_map !$pragma C (r_map)
	external get_str !$pragma C (get_str)

	write(*,*)'----Input Mono-Image File Name---->'
              read(*,*)fill
              call r_map(fill,tmp)

              do i=1,512
              do j=1,512
                 gy(j,i)=tmp(j,513-i)
              end do
              end do
              do i=1,50
                 boxtext(i)=' '
              end do  

              nbox=0
  100         write(*,*)'DO YOU WANT TO MARK LEBALS ON MAP?<Y/N>'
              call get_str(key1,textlen)
             if(key1.eq.'y'.or.key1.eq.'Y') then
	         write(*,*)'You want to Select a Clip Area?<Y/N>'
		 read(*,*)key
		 if(key.eq.'y'.or.key.eq.'Y') then
	            write(*,*)'Please input (I1,I2,J1,J2)[0-512]:'
		    read(*,*)I1,I2,J1,J2
		 endif 
                 WHITE=255.
                 BLACK=0.    
                 call pgbegin(0,'?',1,1)
                 call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
                 call pgbox('bc',0.0,0,'bc',0.0,0)
                 call pgwindow(1.,512.,1.,512.)
	         call gwhbc(gy,tmp,IDIM,JDIM,WHITE,BLACK)
	         CALL GGRAY(tmp,IDIM,JDIM,I1,I2,J1,J2,WHITE,BLACK,TR,box,nbox)

                 write(*,*)'How many strings do you expect to lebal?' 
                 read(*,*)nbox 
  200            write(*,*)'label the  string?<y/N>'
	         call get_str(key2,textlen)
	         if(key2.eq.'y') then
	            write(*,10)nbox 
  10                format('Select string No.[1--',I2,']:')
	            read(*,*)ns
	            xbox(ns)=0.0
                    ybox(ns)=0.0
                    nch=0 
	            write(*,*)'Please locate the position you want to mark!'
	            call pgsci(2)
	            call pgolin(1,nch,xbox(ns),ybox(ns),19)
	            call pgtext(xbox(ns)-4.,ybox(ns)-6.,char(48+ns))
                    write(*,*)'Please input the TEXT in box',ns
                    call get_str(boxtext(ns),textlen)
                    boxlen(ns)=textlen(1)
                    n=boxlen(ns)+1
                    boxtext(ns)(n:n)=' '
	            write(*,*)'how many simbols do you input?'
	            read(*,*)simlen(ns)
	            goto 200
	          end if
                  goto 100
	      end if 
	     do k=1,nbox
	        box(1,k)=xbox(k)-8.
	        box(2,k)=box(1,k)+8*simlen(k)
	        box(3,k)=ybox(k)-10.
	        box(4,k)=box(3,k)+22.
	     end do
	     call pgend

	      
              call pgbegin(0,'?',1,1)
              call pgvsize(1.,1.4*5.28+1.,0.3,1.4*3.627+0.3)
C  ----view field of active region(5'.23 X 3'.626)
C------ 35".0/2.08cm   box field(18.70cm, 12.92cm)
              call pgbox('bc',0.0,0,'bc',0.0,0)
              call pgwindow(1.,512.,1.,512.)
c           boxmeg(1)='\\(0274)_S\\u1\\d'
c           boxmeg(2)='S\\u1\\d..\\d1\\u'
	      call pgslw(1)
	      CALL GGRAY(tmp,IDIM,JDIM,I1,I2,J1,J2,WHITE,BLACK,TR,box,nbox)
              call pgscf(2)
              call pgsch(0.8)
              call pgslw(2)
            do k=1,nbox
             write(*,*)boxtext(k),"length==>",textlen(1)
             call pgtext(box(1,k)+2.,box(3,k)+8.5,boxtext(k))
            end do
              meg='Oct. 28, 1994 \\A \\(0252)  S\\u1\\d   N\\d2\\u'
c             call pgtext(30.,530.,meg)
	      call pgend 
 
	      write(*,*)'quit or redisplay?<q/r>'
	      call get_str(key3,textlen)
	      if(key3.eq.'r') then
	        goto 100
	      end if
	      END

C-------FSub.17
C------ Select an Area with Cursor for Calculating Current Intensity
	SUBROUTINE crnt(density,current,jx1,jx2,iy,num,color,vmin)
	real density(170,128),current,vmin
	integer color
        INTEGER jx1(200),jx2(200),iy(200)
 
 	curt=0.0
	dx=5.23*60./170.
	dy=3.626*60./128.

        do m=1,num
          i=iy(m)
          jj1=jx1(m)
          jj2=jx2(m)
          do j=jj1,jj2
	   if(color.eq.7) then
	     if(density(j,i).gt.vmin) then
	       curt=curt+density(j,i)*1.e-5
 	       yi=float(i)
	       xj=float(j)
	       call pgsci(color)
	       call pgpoint(1,xj,yi,17)
	     end if
	   else if(color.eq.2) then
             if(density(j,i).lt.-vmin) then
               curt=curt+density(j,i)*1.e-5
               yi=float(i)
               xj=float(j)
               call pgsci(color)
               call pgpoint(1,xj,yi,17)
             end if
	   end if
          end do
        end do
	current=curt*dx*dy*(7.26e5)**2
c	write(*,11)current
c11     format('current---> ',e11.4,' (Amp)')
        RETURN
	END
