C------------- FSub. 1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC                                                                           CC
CC               THE   PROGRAM FOR  VERTICAL CURRENT                         CC
CC                                                                           CC
CC                            Wang Tong Jiang                                CC
CC                                                                           CC
CC                               ---- Jun 20,1993 (From Nanjin University)   CC
CC				-----Jan 25,1996			     CC
CC                                                                           CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	Subroutine sub_crt 

        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 bxr(128,128),bxi(128,128),byr(128,128),byi(128,128)
	real jzr(128,128),jzi(128,128)
	real disp(170,128),flz(170,128),rdat(170,128),adat(170,128)
	integer jx1(200),jx2(200),iy(200)
	character file*12,filx*12,fily*12,key,kyy,ckey
	character crm*4,str*60,date*20,cht(8)*5
	character fill*12,fils*12
	integer anw
	character*6 xx,yy,bbc,l,al,bt,fr
	character see*2,time*8
	CHARACTER MESG*400
	real posit(7)
	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),mvdx,mvdy
	real gm(6),trmv(6)
	real trsp(6)

	data gm/90.,140.,170.,200.,240.,256./
	data cm/200.,400.,800.,1000.,1200.,1400.,1600.,2000./
	data cht/'0.20','0.40','0.80','1.00','1.20','1.40','1.60','2.00'/
	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 trsp/0.0,1.0,0.0,0.0,0.0,1.0/	

	external r_map !$pragma C (r_map)
	external r_head !$pragma C (r_head)
	external rw_fit !$pragma C (rw_fit)
	external exel !$pragma C (exel)
	pi=3.1415926

	do m=1,8
	   cp(m)=-cm(m)
	   fp(m)=-fm(m)
	end do

111	write(*,*)'OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO'
	write(*,*)'OOOO      Menu 2 FOR CALCULATION OF CURRENT             OOOO'
	write(*,*)'OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO'
	write(*,*)'OOOO   0. goto ----> <Main Menu>                        OOOO'
	write(*,*)'OOOO   1. read bi*.dat in fits format for current data  OOOO'
	write(*,*)'OOOO   2. the difference method for current             OOOO'
	write(*,*)'OOOO   3. the fourier method for current                OOOO'
	write(*,*)'OOOO   4. display or print for current                  OOOO'
	write(*,*)'OOOO   5. calculate current intensity for an area       OOOO'
	write(*,*)'OOOO   6. lap contour-map on gray-map                   OOOO'
	write(*,*)'OOOO   7. make magnetic-map location same as monoimage  OOOO'
	write(*,*)'OOOO   8. quit                                          OOOO'
	write(*,*)'OOOO   9.<Help>                                         OOOO'
	write(*,*)'OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO'
	
	write(*,*)'######## SELECT NUMBER:(0-9)===>  '
	READ(*,*)ANW 

	IF(ANW.EQ.0) THEN
	   RETURN
	ELSE IF(ANW.EQ.1) THEN
	   write(*,*)'PLEASE INPUT bi*.dat====>'
	   read(*,*)file
	   call rw_fit(file,disp,mesg,'r')
	   rmax=-1.e6
	   rmin=1.e6
	   do i=1,128
	   do j=1,170
	      if(rmax.lt.disp(j,i)) rmax=disp(j,i)
	      if(rmin.gt.disp(j,i)) rmin=disp(j,i)
           end do
	   end do
	   write(*,5)rmax*1.e-5,rmin*1.e-5

	ELSE IF(ANW.EQ.2.OR.ANW.EQ.3) THEN
	   write(*,*)'PLEASE ENTER INPUT BX-FILENAME===>'
	   read(*,*)filx
	   call rw_fit(filx,bx,mesg,'r')
           call FMAXMIN(bx,170,128,bxmax,bxmin)
           write(*,*)'bxmax==> bxmin==>',bxmax,bxmin
	   write(*,*)'PLEASE ENTER INPUT BY-FILENAME===>'
	   read(*,*)fily
	   call rw_fit(fily,by,mesg,'r')
	   call FMAXMIN(by,170,128,bymax,bymin)
	   write(*,*)'bymax==> bymin==>',bymax,bymin

C   ------- BEGIN TO CALCULATE THE CURRENT -------
	  dx=5.23*60./170.
	  dy=3.626*60./128.
C-------- dx and dy are angle-second per pix in x-axis and in y-axis

C----------the difference method for current--------
	   if(anw.eq.2) then
	     write(*,*)'------the difference method---BEGIN......'
	     cail=1.e-4/(7.26e5*4.*pi*1.e-7)
	     rmax=-1.e6
	     rmin=1.e6
	     do i=2,127
	     do j=2,169
	        dxx=(bx(j,i+1)-bx(j,i-1))/(2.0*dy)
	        dyy=(by(j+1,i)-by(j-1,i))/(2.0*dx) 
	        disp(j,i)=(dyy-dxx)*cail*1.e5
C          -------- Cailbration about current ----------
	        if(rmax.lt.disp(j,i)) rmax=disp(j,i)
	        if(rmin.gt.disp(j,i)) rmin=disp(j,i)
             end do
	     end do
	     write(*,5)rmax*1.e-5,rmin*1.e-5
5	     format('current---max: ',e11.4,'A/(m*m) ',e11.4,'A/(m*m)')

C      --------the fourier method for current -----------
	  else if(anw.eq.3) then
	    write(*,*)'-------the fourier method---BEGIN.......'
	    do i=1,128
	    do j=1,128
	       jl=j+20
	       bxr(j,i)=bx(jl,i)
	       bxi(j,i)=0.0
	       byr(j,i)=by(jl,i)
	       byi(j,i)=0.0
	    end do
	    end do	  
	    call d2ft(7,128,1.,bxr,bxi)
	    call d2ft(7,128,1.,byr,byi)
	    do m=1,128
	    do n=1,128
	     jzr(m,n)=-byi(m,n)*sin(2*pi*m/128.)/dx+bxi(m,n)*sin(2*pi*n/128.)/dy
	     jzi(m,n)=byr(m,n)*sin(2*pi*m/128.)/dx-bxr(m,n)*sin(2*pi*n/128.)/dy
	    end do
	    end do
	    call d2ft(7,128,-1.,jzr,jzi)
	    cail=1.e-4/(7.26e5*4.*pi*1.e-7)
	    do i=1,128
	    do j=1,128
	    jzr(j,i)=1.e5*cail*jzr(j,i)
	    end do
	    end do
	     do i=10,120
	     do j=10,120
	        jl=j+20
	        jzr(j,i)=-jzr(j,i)
	        disp(jl,i)=jzr(j,i)
	        if(rmax.lt.jzr(j,i)) rmax=jzr(j,i)
	        if(rmin.gt.jzr(j,i)) rmin=jzr(j,i)
	        if(zima.lt.jzi(j,i)) zima=jzi(j,i)
	        if(zimi.gt.jzi(j,i)) zimi=jzi(j,i)
	     end do
	     end do
	   end if

             write(*,*)'--- Output File Name ------'
             read(*,*)file
             call rw_fit(file,disp,mesg,'w')

C------- Dispaly Current with Contour Map or Gray Map
	ELSE IF(ANW.EQ.4.or.ANW.EQ.5) THEN
	   write(*,*)'bimax==> bimin==> ',rmax,rmin
 	   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)

 	   write(*,*)'Draw Gray-map or Contour-map? <G/C>'
	   read(*,*)key
	   if(key.eq.'g'.or.key.eq.'G') then
	     rat=0.7
	     write(*,*)'Change Gray-level [0.7]?<Y/N>'
	     read(*,*)key
	     if(key.eq.'y') then
	       write(*,*)'Please input New level (0.0--1.0):'
	       read(*,*)rat
	     end if
             call pggray(disp,170,128,1,170,1,128,rmin*rat,rmax*rat,tr)
           else if(key.eq.'c'.or.key.eq.'C') then
             CALL GSLS(2)
             call pgsci(2)
             CALL PGSLW(1)
             CALL PGCONT(disp,170,128,1,170,1,128,CP,-8,TR)
             call pgsls(1)
             call pgsci(7)
             CALL PGSLW(1)
             CALL PGCONT(disp,170,128,1,170,1,128,CM,-8,TR)
	   end if

C--------- Calculate Current Intensity for an Area
	   if(anw.eq.5) then
	       Current=0.0
	       Ptotal=0.0
	       Gtotal=0.0
888	       write(*,*)'--Give the Min Limit of Integration ==>[400]:'
	       write(*,*)'-- (value 400 , i.e. 0.004 Am/m*m )' 
               read(*,*)vmin
	       WRITE(*,*)'Calculate Positive-I(Up) or Negative-I(Down) ?<P/N>'
	       read(*,*)key
	       if(key.eq.'p') then
	      	  call boundary(jx1,jx2,iy,num)
	          call crnt(disp,current,jx1,jx2,iy,num,7,vmin)
	          ptotal=ptotal+current
		  write(*,25)current
25		  format('Positive current==>',e11.4,'Amp')
	       else if(key.eq.'n') then
		  call boundary(jx1,jx2,iy,num)
		  call crnt(disp,current,jx1,jx2,iy,num,2,vmin)
	  	  gtotal=gtotal+current
		  write(*,15)current
15	          format('Negative current==>',e11.4,'Amp')
	       end if
	       write(*,35)ptotal
35	       format('Sum (+)Current:',e11.4,'(Amp)')
	       write(*,45)gtotal
45	       format('Sum (-)Current:',e11.4,'(Amp)')
	       write(*,*)'Caculate Current for next Area?<Y/N>'
	       read(*,*)key
	       if(key.eq.'y') goto 888
	   end if
	   call pgend  
	   
C------- Lap Contour Map on Gray Map
	ELSE IF(ANW.EQ.6) THEN
	   write(*,*)'Select Gray-map Type: MONO-IMAGE or CURRENT-IMAGE?<M/C>'
	   read(*,*)key
	   if(key.eq.'m') then
	      write(*,*)'----Input Mono-Image File Name---->'
              read(*,*)fill
              call r_map(fill,tmp)
 	      call r_head(fill,date,time,see,al,bt,bbc,l,xx,yy,fr,posit)
	      date=time//' '//date
	      write(*,*)'date== ',date

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

	      WHITE=255.
	      BLACK=0.
	      WRITE(*,*)'--- Change the Contrast of Gray?<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,512
	         do j=1,512
	            gy(j,i)=(gy(j,i)+gy0)/255.
	            gymin=min(gy(j,i),gymin)
	            gymax=max(gy(j,i),gymax)
	         end do
	         end do
	         WHITE=gymax*ratio
	         BLACK=gymin*ratio
	      END IF

	      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,TR)

              call pgwnad(0.,105.,0.,78.)
              call pgvsize(0.05,10.5+0.05,0.05,7.8+0.05)
              call pgscf(2)
	      call pgsch(0.8)
	      call pgslw(1)
              str=fill//' Time: '//date
              call pgtext(15.,61.,str)
              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)

	   else if(key.eq.'c') then
	      rat=1.0
	      write(*,*)'--- Input Black-White Contrast of Gray-Map==>[1.0]'
	      read(*,*)rat
	      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.,170.,1.,128.)
	      call pggray(disp,170,128,1,170,1,128,rmin*rat,rmax*rat,tr)
              call pgwnad(0.,105.,0.,78.)
              call pgvsize(0.05,10.5+0.05,0.05,7.8+0.05)
              call pgscf(2)
              call pgsch(0.8)
              str='the Vertical Current Gray-Map '
              call pgtext(15.,65.,str)
              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)
	   end if

	   write(*,*)'Select Coutour-Map type: Longitude-field or Current?<L/C>'
	   read(*,*)key

	   write(*,*)'Change Position of Second Map?<Y/N>'
	   read(*,*)kyy
	   if(kyy.eq.'y') then
	     tr(1)=trmv(1)
	     tr(4)=trmv(4) 
	   else
	     tr(1)=0.0
	     tr(4)=0.0
	   end if

             write(*,*)'You hope the White line on gray image<Y/N>?'
             read(*,*)ckey

	   if(key.eq.'l') then
	     write(*,*)'-----Input Longitude Field File:[*.bzo]--->'
	     read(*,*)fils
	     XXG=GTOTAL
	     call rw_fit(fils,flz,mesg,'r')
	     GTOTAL=XXG
C------------COPY GTOTAL TO XXG BEFORE CALL RW_FIT,BECAUSE IT DAMAGE GTOTAL
             call pgwindow(1.,170.,1.,128.)
	     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(1)
	        call pgsci(10)
	        if(ckey.eq.'y'.or.ckey.eq.'Y') then
		   call pgsci(0)
	  	endif
		CALL MVARROW(RDAT,ADAT,1,170,1,128,dqu,tr(1),tr(4))
	    END IF

	     
	     CALL GSLS(2)
             call pgsci(2)
             if(ckey.eq.'y'.or.ckey.eq.'Y') then 
                   call pgsci(0) 
             endif
             CALL PGSLW(1)
             CALL PGCONT(flz,170,128,1,170,1,128,FP,-8,TR)
             call pgsls(1)
             call pgsci(7)
	     call pgslw(3)
             if(ckey.eq.'y'.or.ckey.eq.'Y') then 
                   call pgsci(0) 
		   call pgslw(1)
             endif
             CALL PGSLW(1)
             CALL PGCONT(flz,170,128,1,170,1,128,FM,-8,TR)

             call pgwnad(0.,105.,0.,78.)
             call pgvsize(0.05,10.5+0.05,0.05,7.8+0.05)
             call pgscf(2)
	     call pgslw(1)
	     call pgsch(0.8)
	     str='Longitude Magnetic field: '//fils
	     call pgtext(15.,63.,str)
             call pgtext(24.,59.,'2000 Guass')
             call pgtext(42.,59.,'1000 Guass')
             call pgend
       WRITE(*,*)'GTOTAL=>',GTOTAL
       WRITE(*,*)'PTOTAL==>',PTOTAL

 	   else if(key.eq.'c') then
	     call pgwindow(1.,170.,1.,128.)
             CALL GSLS(2)
             call pgsci(2)
             if(ckey.eq.'y'.or.ckey.eq.'Y') then
                   call pgsci(0)
             endif
             CALL PGSLW(1)
             CALL PGCONT(disp,170,128,1,170,1,128,CP,-8,TR)
             call pgsls(1)
             call pgsci(7)
	     call pgslw(3)
             if(ckey.eq.'y'.or.ckey.eq.'Y') then
                   call pgsci(0)
	           call pgslw(1)
             endif
             CALL PGCONT(disp,170,128,1,170,1,128,CM,-8,TR)

             call pgwnad(0.,105.,0.,78.)
             call pgvsize(0.05,10.5+0.05,0.05,7.8+0.05)
             call pgscf(2)
             call pgsch(2.)
	     call pgslw(1)
             call pgtext(10.,73.,'Vertical Current of Solor Disk ')
             call pgsch(0.8)
	     call string(rmax,crm)
	     str='Up    Current Density Max: '//crm(1:1)//'.'//crm(2:4)//
     1	         '\\(0727)10\\u-2\\d (A/m\\u2\\d)'
	     call pgtext(15.,69.,str)
	     call string(-rmin,crm)
	     str='Down Current Density Max: '//crm(1:1)//'.'//crm(2:4)// 
     1           '\\(0727)10\\u-2\\d (A/m\\u2\\d)'
	     call pgtext(15.,67.,str)
	     pto=ptotal*1.e-10
	     call string(pto,crm)
	     if(crm(1:1).eq.'0') crm(1:1)=' '
	     str='Up    Current Intensity:  '//crm(1:2)//'.'//crm(3:4)// 
     1           '\\(0727)10\\u12\\d (A)'
	     call pgtext(15.,65.,str)
	     gto=gtotal*1.e-10
	     call string(-gto,crm)
             if(crm(1:1).eq.'0') crm(1:1)=' ' 
	     str='Down Current Intensity:  '//crm(1:2)//'.'//crm(3:4)// 
     1           '\\(0727)10\\u12\\d (A)'
	     call pgtext(15.,63.,str)
c	     str=fill//' Time: '//date
c 	     call pgtext(15.,61.,str)
	     call string(vmin,crm)
	     str='Integration for Current'
	     call pgtext(60.,69.,str)
	     str='the Min Limit:'
	     call pgtext(60.,67.,str)
	     str=crm(1:1)//'.'//crm(2:4)//'\\(0727)10\\u-2\\d (A/m\\u2\\d)'
	     call pgtext(60.,65.,str) 

             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-1,yy2-8,'Up   : \\(0796)\\(0796)\\(0796)\\(0796)')
             call pgtext(xx2-1,yy2-10.,'Down : ------')
             call pgtext(xx2-1,yy2-12.,'Level : A/m\\u2\\d')
	     call pgtext(xx2-1,yy2-14.,'(\\(0727)10\\u-2\\d)')
             do i=1,8
                call pgtext(xx2+2.,yy2-14.-2.6*i,cht(i))
             end do
	     call pgend
	   end if
	   tr(1)=0.0
	   tr(4)=0.0


	ELSE IF(ANW.EQ.7) THEN
	      write(*,*)'----Input Mono-Image File Name---->'
              read(*,*)file
              call r_map(file,tmp)
              do i=1,512
              do j=1,512
                 gy(j,i)=tmp(j,513-i)
              end do
              end do

             write(*,*)'-----Input Longitude Field File:[*.bzo]--->'
             read(*,*)file
             call rw_fit(file,flz,mesg,'r')
	     trmv(1)=0.0
	     trmv(4)=0.0

777	      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 pgwindow(1.,512.,1.,512.)
             call pgsls(1)
             call pgsci(5)
             CALL PGSLW(2)
             CALL PGCONT(gy,512,512,1,512,1,512,GM,-6,TR)

             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(flz,170,128,1,170,1,128,FP,-6,TRmv)
             call pgsls(1)
             call pgsci(7)
             CALL PGSLW(1)
             CALL PGCONT(flz,170,128,1,170,1,128,FM,-6,TRmv)
	     write(*,*)'Change Position of Longitude-field Map?<Y/N>'
	     read(*,*)key
	     if(key.eq.'y') then 
		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  mono-image map (BLUE)--'
	        call pgsci(5)
                call pgolin(20,nch,xmo,ymo,18)
	 	   
	write(*,*)'---Select position  points on longitude-field (YELLOW-RED)--'
	        nch=0
	        call pgsci(4)
	        call pgolin(20,nch,xmg,ymg,17)
	        write(*,*)'Select point number as contrast:-->',nch
	        mvdx=0.0
	        mvdy=0.0
	        do k=1,nch
	           mvdx=mvdx+(xmo(k)-xmg(k))
		   mvdy=mvdy+(ymo(k)-ymg(k))
	        end do
	        mvdx=mvdx/nch
	        mvdy=mvdy/nch
	        write(*,55)mvdx,mvdy
55  	        format('move dx= ',f8.3,'  move dy= ',f8.3)
	        trmv(1)=trmv(1)+mvdx
	        trmv(4)=trmv(4)+mvdy
		call pgend
		goto 777
	    else if(key.eq.'n') then
	      call pgend
	      write(*,*)'you Affirm the result or Restore origin pattern?<A/R>'
	      read(*,*)kyy
		if(kyy.eq.'r') then
		  trmv(1)=0.0
		  trmv(4)=0.0
	          goto 777
	        end if
	    end if
	      
	ELSE IF(ANW.EQ.8) THEN
	    goto 999
	ELSE IF(ANW.EQ.9) THEN
	    call exel
        END IF
	GO TO 111
999	write(*,*)'------ Bye !  Thank You for Using ! ------'
        stop
        end


