CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC                                                                          CC
CC     Judge the Direction of the Transversional Magnetic Field             CC
CC                                                                          CC
CC                         WANG TONG JIANG                                  CC
CC                                                                          CC
CC                              ----- Jun 5,1993 (From Nanjin University)   CC 
CC             ****** Version 1.1**** (MAY 6, 1994 in Nanjin University)    CC
CC             --     Version 1.2  --  (Jan 25, 1996 Huairou Station )      CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC                                                                          CC
CC     THIS PROGRAM IS ABOUT THE JUDGEMENT OF THE 180 DEGREE AMBIGUITY  OF  CC
CC  THE TRANSVERTIONAL MAGNETIC FIELD.THE PARTS OF SUBROUTINE IS PROVIDED   CC
CC  BY PROFESSOR WANG JINXIOU WHO DO WELL IN THE ASPECT,AND GAVE ME MORE    CC
CC  GIUDENCE . THE CONTENT IN DIFFERENCE OF HIS PROGRAM  IS THAT FIRST,THE  CC
CC  SMOOTH OF 3 X 4 IS DONE BEFORE THE JUDGEMENT AND SECOND,THE MAGNTIC     CC
CC  LINE METHOD AND THE BOX CURSOR METHOD ARE ADDED UP. IN ADDITION,TWO     CC
CC  KINDS OF THE CAILBRATION METHODS OF TRANSVERSE FIELD AND LONGITUDE      CC
CC  FIELD  ARE REFERRED. THE OBTAINED RESULTS OF VECTOR MAGNETIC FIELD ARE  CC
CC  SAVED IN THREE DIFFERENT FILES WITH THE FITS FORMAT. THE UNIT OF DATA   CC
CC  IN THREE FILES IS GAUSS.                                                CC
CC                                                                          CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        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 gdat(170,128),RDAT(170,128),adat(170,128),TR(6),CM(8),CP(8)
	character bznam*13,temp*9
	real padat(170,128),prdat(170,128)
	real nbx(150),nby(150),bb(150)
	integer anw
	character key,CHT(8)*5,bznam*13,bynam*13,bxnam*13
	CHARACTER CDATE*30,SEEING*15,WAVE*25,CFILE*51,CONT*35
        CHARACTER CARR*40,FRAME*30,SUNDISK*40,ORTHO*40
	character crm*4,str*52
	real xdat(170,128),ydat(170,128),zdat(170,128)
        integer jx1(200),jx2(200),iy(200)
	
        CHARACTER MESG*400

        DATA TR/0.0,1.0,0.0,0.0,0.0,1.0/
	data cm/-2880.0,-2560.0,-2240.,-1920.,-1280.,-640.,-160.,-40./
	DATA CHT/'2880.','2560.','2240.','1920.','1280.','640.','160.','40.'/

	external rw_fit !$pragma C (rw_fit)
	external exel !$pragma C (exel)
	
	WRITE(*,*)' '
	WRITE(*,*)' ________________________________________________________'
	WRITE(*,*)'- Software of Solar Vector Magnetic Field Process     ---'
	WRITE(*,*)'-            (FeI5324.19A)                            ---'
	WRITE(*,*)'-      Designer   Wang Tongjiang	                 ---'
	WRITE(*,*)'-      Instructor Ai Guoxiang                         ---'
	WRITE(*,*)'-          Version 1.2 (1993--1996)                   ---'
	WRITE(*,*)'_________________________________________________________'

C------ contour level of longitudial magnetic field
	do k=1,8
	   cp(k)=-cm(k)
	end do
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--- JUDGEMENT ABOUT THE 180-DEGREE AMBIGOUSE OF THE TRANSVERTIONAL FIELD----- 

100	write(*,*)' '
	write(*,*)' '
	write(*,*)'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
	write(*,*)'WWWWW        < THE MAIN MENU >                         WWWW'
	write(*,*)'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
	write(*,*)'WW  0.  Calibration for Stokes L, Q and U              WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  1.  display the observed vector magnetic field     WWWW'
	write(*,*)'WW  2.  display vector field with direction            WWWW'
	write(*,*)'WW  3.  display the calculated potential field         WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  4.  adjust local area  by cursor                   WWWW'
	write(*,*)'WW  5.  adjust box-area by cursor                      WWWW'
	write(*,*)'WW  6.  modify field direction by magnetic line        WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  7.  read data(*.bxo,*.byo,*.bzo)  and display      WWWW'
	write(*,*)'WW  8.  save the Observed Vector Field                 WWWW'
	write(*,*)'WW  9.  save the potential field                       WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  11. [Menu 1: Direction]                            WWWW'
	write(*,*)'WW  12. [Menu 2: Current  ]                            WWWW'
	write(*,*)'WW  13. [Menu 3: Project  ]                            WWWW'
	write(*,*)'WW  14. [Menu 4: Spectrum ]                            WWWW'
	write(*,*)'WW  15. [Menu 5: Function ]                            WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  10. Quit                                           WWWW'
	write(*,*)'WW --------------------------------------------------- WWWW'
	write(*,*)'WW  20. Help                                           WWWW'
	write(*,*)'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
	write(*,*)'WW      PLEASE ENTER THE SELECTED NUMBER!              WWWW'
	write(*,*)'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
	write(*,*)' '
	write(*,*)'###### NO.(0-16)===>  '
	read(*,*)anw
	if(anw.eq.1.or.anw.eq.2.or.anw.eq.3.or.anw.eq.4.or.anw.eq.5.or.
     1      anw.eq.6.or.anw.eq.7) then
		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)
C------- the new data come from observation June 1995 by Zhang H.
                call pgwindow(x1,x2,y1,y2)
                call pgbox('bcs',0.0,0,'bcs',0.0,0)

                xlen=169.*35./(5.26*60.)
                call pgsch(0.7)
                call pgptext(1.0,133.,0.0,0.0,'|')
                call pgptext(1.0+xlen,133.,0.0,0.0,'|')
                dxlen=xlen/20.
                do i=0,18
                   zxx=1.0+dxlen*float(i)
                   call  pgptext(zxx,133.,0.0,0.0,'-')
                enddo
                call pgscf(2)
                call pgptext(xlen/2.,136.,0.0,0.5,'35\\u"\\d.0')
                call pgsch(1.)
	 endif

	    IF(ANW.EQ.0) THEN
	        CALL calib(gdat,rdat,adat,dqu,mesg)
	        cdate=MESG(30:60)
	        cfile=MESG(60:110)
                sundisk=MESG(110:150)
                ortho=MESG(150:190)
                carr=MESG(190:230)
                frame=MESG(230:260)
                seeing=MESG(260:280)
		cont=MESG(280:320)
		wave=MESG(320:350)
	        temp=cfile(13:21)
		goto 100
	     ELSE IF(ANW.EQ.1) THEN
 		CALL obs_arr(rdat,adat,nx1,nx2,ny1,ny2,dqu)
	     ELSE IF(ANW.EQ.2) THEN
	        CALL  arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
	     ELSE IF(ANW.EQ.3) THEN
		CALL  arrow(prdat,padat,nx1,nx2,ny1,ny2,tmax)
	        write(*,*)'potential field Btmax-->',tmax
	     ELSE IF(ANW.EQ.4) THEN
                call arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
  	        call boundary(jx1,jx2,iy,num)
		CALL adcur(rdat,adat,jx1,jx2,iy,num)
	     ELSE  IF(ANW.EQ.5.OR.ANW.EQ.6) THEN
C------------------------ Draw the Grid in Frame -----------------------
		  DO I=0,130,10
		  DO J=1,170
		     gx=float(j)
		     gy=float(i)
		     call pgpoint(1,gx,gy,-1)
		  END DO 
		  END DO
		  DO I=1,128
		  DO J=0,170,10
		     gx=float(j) 
                     gy=float(i) 
                     call pgpoint(1,gx,gy,-1) 
                  END DO  
                  END DO 
                  IF(ANW.EQ.5) THEN
                     call arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
                     call box_adcur(rdat,adat,nbx,nby,bb,ll)
                     call boxarr(bb,nbx,nby,ll)
                  ELSE IF(ANW.EQ.6) THEN
                     CALL new_magline(rdat,adat,padat,bb,nbx,nby,ll)
                     call arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
                     call boxarr(bb,nbx,nby,ll)
	          END IF
	     ELSE IF(ANW.EQ.7) THEN
		   write(*,*)'input file bx-name==>'
		   read(*,*)bxnam
                   call rw_fit(bxnam,xdat,MESG,'r')
		   write(*,*)'input file by-name==>'
		   read(*,*)bynam
                   call rw_fit(bynam,ydat,MESG,'r')
		   write(*,*)'input file bz-name==>'
		   read(*,*)bznam
                   call rw_fit(bznam,gdat,MESG,'r')
	           cdate=MESG(30:60)
	           cfile=MESG(60:110)
                   sundisk=MESG(110:150)
                   ortho=MESG(150:190)
                   carr=MESG(190:230)
                   frame=MESG(230:260)
                   seeing=MESG(260:280)
		   cont=MESG(280:320)
		   wave=MESG(320:350)
	           temp=bznam(1:9)
		   rmax=0.0
		   gmax=-1.e6
		   gmin=1.e6
		   do i=1,128
	 	   do j=1,170
	 	      rdat(j,i)=sqrt(xdat(j,i)**2+ydat(j,i)**2)
		      adat(j,i)=atan2d(ydat(j,i),xdat(j,i))
		      if(rmax.lt.rdat(j,i)) rmax=rdat(j,i)
           	      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)')
		   write(*,10)rmax
10                 format('btmax==>',f8.2,' Guass')
		   dqu=1.2*rmax
		   call arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
	     ELSE IF(ANW.EQ.8) THEN
                   do i=1,128
                   do j=1,170
                       xdat(j,i)=rdat(j,i)*cosd(adat(j,i))
                       ydat(j,i)=rdat(j,i)*sind(adat(j,i))
		       zdat(j,i)=gdat(j,i)
                   end do
                   end do
                   bxnam=temp(1:9)//'bxo '
                   bynam=temp(1:9)//'byo '
                   bznam=temp(1:9)//'bzo '
C -----   DATA TYPE IN OUTPUT FILE: SHORT ------
		   call rw_fit(bznam,zdat,MESG,'w')
	           call rw_fit(bxnam,xdat,MESG,'w')
		   call rw_fit(bynam,ydat,MESG,'w')
	           goto 100
	    ELSE IF(ANW.EQ.9) THEN
                   do i=1,128
                   do j=1,170
                       xdat(j,i)=prdat(j,i)*cosd(padat(j,i))
                       ydat(j,i)=prdat(j,i)*sind(padat(j,i))
		       zdat(j,i)=gdat(j,i)
                   end do
                   end do
                   bxnam=temp(1:9)//'bxp '
                   bynam=temp(1:9)//'byp '
                   bznam=temp(1:9)//'bzp '
C -----   DATA TYPE IN OUTPUT FILE: SHORT ------
		   call rw_fit(bznam,zdat,MESG,'w')
	           call rw_fit(bxnam,xdat,MESG,'w')
		   call rw_fit(bynam,ydat,MESG,'w')
	           goto 100
	    ELSE IF(ANW.EQ.10) THEN
		   GOTO 999
	    ELSE IF(ANW.EQ.20) THEN
	        call exel    ! C routine for more help
	        goto 100
            ELSE IF(ANW.EQ.11) THEN
	        goto 111   !--- Direction
	    ELSE IF(ANW.EQ.15) THEN
		goto 222   !--- Function 
	    ELSE IF(ANW.EQ.12) THEN
	        call sub_crt  
	        goto 100
	    ELSE IF(ANW.EQ.13) THEN
		call sub_prj 
	        goto 100
	    ELSE IF(ANW.EQ.14) THEN
		call spec_mag
	        goto 100
	    END IF

            CALL GSLS(2)
	    call pgsci(2)
            CALL PGSLW(1)
            CALL PGCONT(gdat,170,128,1,170,1,128,CM,-8,TR)
            call pgsls(1)
	    call pgsci(7)
            CALL PGSLW(1)
            CALL PGCONT(gdat,170,128,1,170,1,128,CP,-8,TR)
        call pgsci(7)
        call pgwnad(0.,105.,0.,100.)
        call pgvsize(0.05,10.5+0.05,0.05,10.0+0.05)
        call pgscf(2)
        call pgsch(2.)
         call pgtext(10.,73.,'SOLOR VECTOR MAGNETIC FILED ')
        call pgsch(0.8)
        call pgtext(15.,69.,cdate)
        call pgtext(15.,67.,cfile)
c       call pgtext(15.,65.,cont)
        call pgtext(15.,63.,sundisk)
c       call pgtext(15.,61.,ortho)
        call pgtext(15.,65.,Carr)
        call pgtext(72.,69.,wave)
        call pgtext(72.,67.,seeing)
        call pgtext(72.,65.,Frame)
        call pgline(3,nxx,nyy)
        call pgtext(xx1-0.7705,yy1,'\\(2247)')
        call pgtext(xx1-0.775,yy1+2.,'N')
        call pgtext(xx3,yy3-.5,'> W')
        call pgtext(xx2,yy2-8,'N : \\(0796)\\(0796)\\(0796)\\(0796)')
        call pgtext(xx2,yy2-10.,'S : ------')
        call pgtext(xx2,yy2-12.,'Lev. : Gauss')

	CALL STRING(GMAX,CRM)
	STR='Bz (+N) Max: '//crm//' Guass'
	CALL STRING(-GMIN,CRM)
	STR=STR(1:27)//' Bz (-S) Max: '//crm//' Guass'
	call pgtext(15.,61.,str)
	CALL STRING(RMAX,CRM)
	STR='Bt-Field Max: '//crm//' Guass'
	call pgtext(15.,59.,str)
	call pgtext(24.,55.,'2000 Guass')
	call pgtext(42.,55.,'1000 Guass')
        do i=1,8
          call pgtext(xx2+2.,yy2-14.-2.6*i,cht(i))
        end do
	call pgiden
	call pgend
	GOTO 100

	
111	write(*,*)' '
	write(*,*)' '
	write(*,*)' '
	write(*,*)'*******************************************************'
	write(*,*)'********** < Menu 1 for Judging Direction > ***********'
	write(*,*)'*******************************************************'
	write(*,*)'*   0.  GO ON .......return to Main Menu              *'
	write(*,*)'*   1.  magnetic charge method                        *'
	write(*,*)'*   2.  potential field method                        *'
	write(*,*)'*   3.  krall method                                  *'
	write(*,*)'*   4.  diverge method                                *'
	write(*,*)'*   5.  concord ratio of Bz/Bt by Potential Field     *'
	write(*,*)'*   6.  Quit                                          *'
	write(*,*)'*******************************************************'
	write(*,*)'*** PLEASE ENTER THE SELECTED NUMBER!******************'
	write(*,*)'*******************************************************'
        write(*,*)'######  NO.(0-6)==>  '
	read(*,*)anw

	if(ANW.EQ.1.OR.ANW.EQ.2.OR.ANW.EQ.3.OR.ANW.EQ.4) THEN
	   do m=1,128
	   do n=1,170
	     prdat(n,m)=0.0
	     padat(n,m)=0.0
	   end do
	   end do
	end if

	IF(ANW.EQ.0) THEN
	  GO TO 100
	ELSE IF(ANW.EQ.1) THEN
	   write(*,*)'   1.  magnetic charge method BEGIN...............'
	   call charge(gdat,rdat,adat,prdat,padat,tmax)
	ELSE IF(ANW.EQ.2) THEN
	   write(*,*)'    2.  potential field method BEGIN...............'
	   call poten(gdat,rdat,adat,prdat,padat,tmax)
	ELSE IF(ANW.EQ.3) THEN
	   write(*,*)'    3.  krall method BEGIN.........................'
	   call krall(gdat,rdat,adat,prdat,padat,tmax)
	ELSE IF(ANW.EQ.4) THEN
	   write(*,*)'   4.  diverge method BEGIN.......................'
	   call diver(gdat,rdat,adat)
	ELSE IF(ANW.EQ.5) THEN
	   rat=rmax/tmax
	   do i=1,128
	   do j=1,170
	      gdat(j,i)=rat*gdat(j,i)
	   end do
	   end do 
	   gmax=rat*gmax
	   gmin=rat*gmin
	ELSE IF(ANW.EQ.6) THEN
	   GOTO 999
	END IF
	   write(*,*)'tmax===>',tmax
  	   write(*,*)'dqu= ',dqu
	GOTO 100

C----------------- Add Menu -----------------------------------------------
222	write(*,*)' '
	write(*,*)' '
	write(*,*)' '
	write(*,*)'    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
	write(*,*)'    $             < Menu 5 Function >                    $'
	write(*,*)'    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
	write(*,*)'    $   0.  RETURN ---> <Main Menu>                      $'
	write(*,*)'    $   1.  Flux of Longitude field for a Area           $'
	write(*,*)'    $   2.  Shear Anger for a Area                       $'
	write(*,*)'    $   3.  Transverse field Flux for a Segnment         $'
	write(*,*)'    $   4.  Mark the Labals in a site                    $'
	write(*,*)'    $   5.  Quit                                         $'
	write(*,*)'    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
        write(*,*)'  '
 	write(*,*)'##### PLEASE ENTER THE SELECTED NUMBER:[No. 0--5]'
	read(*,*)ANW
	if(ANW.NE.0.AND.ANW.NE.5.AND.ANW.NE.4) THEN
            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(1.,170.,1.,128.)
            call pgbox('bcs',0.0,0,'bcs',0.0,0)
            CALL GSLS(2)
            call pgsci(2)
            CALL PGSLW(1)
            CALL PGCONT(gdat,170,128,1,170,1,128,CM,-8,TR)
            call pgsls(1)
            call pgsci(7)
            CALL PGSLW(1)
            CALL PGCONT(gdat,170,128,1,170,1,128,CP,-8,TR)
	end if

	IF(ANW.EQ.0) THEN
	    goto 100
	ELSE IF(ANW.EQ.1) THEN
210	    call boundary(jx1,jx2,iy,num)
            call L_flux(gdat,jx1,jx2,iy,num) 
            write(*,*)'---- You want to select next area?<Y/N>'
            read(*,*)key
            if(key.eq.'y') goto 210 
	ELSE IF(ANW.EQ.2) THEN
	    CALL  arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
220	    call boundary(jx1,jx2,iy,num)
 	    call shear(adat,rdat,padat,jx1,jx2,iy,num)
            write(*,*)'---- You want to select next area?<Y/N>'
            read(*,*)key
            if(key.eq.'y') goto 220
	ELSE IF(ANW.EQ.3) THEN
            CALL  arrow(rdat,adat,nx1,nx2,ny1,ny2,dqu)
 	    call T_flux(rdat,adat)
	ELSE IF(ANW.EQ.4) THEN
	    write(*,*)'You select Magnetic field or Gray-image?<M/G>'
	    read(*,*)key
	    if(key.eq.'m'.or.key.eq.'M') then
 	       call labcn(gdat,rdat,adat,dqu,mesg) 
	    else
	       call labgy
	    endif
	    goto 222
	ELSE IF(ANW.EQ.5) THEN
	    goto 999
	END IF	
	CALL PGEND
	GOTO 222
	          

999	write(*,*)'***********  O.K. Thank You for Using !   **************'
	write(*,*)'------- Vision 1.2  Software of Science ----------------'
	write(*,*)'*********** Wang Tongjiang on 1993-1996 Huairou station***'
	END
	

