
c-----------------------------------------------------------------------
      program  solar
c-----------------------------------------------------------------------
c  traces path of ejected "planet" from surface of sun under gravity,
c    electrostatic and magnetic forces.
c solar2.f same as solar1.f except that uses HAMMING for predictor-corrector 
c    method to control step size and Runge Kutta.
c
c reads physical constants
      common bsun,rsun,rplanet,g,epsilon,sunpot,pie,qsun,qplanet,
     1    rmplanet,rmsun,n,nsteps,ndprint,nprint
      common/out/ sx,sy,velx,vely,gacc,eacc,rmaccx,rmaccy,rmaccz,H,s
      DIMENSION Y(6)
      EXTERNAL F,OUTPUT,FAIL
      integer CFACT
      open(3,file='physconst')
      rewind 3
c radius sun(m),radius planet(m),mass sun(kg), mass planet(kg)
      read(3,*)	
      read(3,*) rsun,rplanet,rmsun,rmplanet
c sun potential(V),bsun is magnetic field at surface of sun at sun's equator
      read(3,*)	
      read(3,*) sunpot,bsun 
c timestep(days),nsteps, velplanet(m/s),start x,y,z(m)
c x,y in plane of equator of sun, ejection initially in x direction
      read(3,*)	
      read(3,*) dt,nsteps,velplanet,sx,sy,sz 
c every nprint steps printed
      read(3,*)	
      read(3,*) ndprint 
      close(3)
c-----------------------------------------------------------------
c Write input values
      open(5,file='record')
      write(*,25) rsun,rplanet,rmsun,rmplanet 
      write(*,26) sunpot,bsun  
      write(5,25) rsun,rplanet,rmsun,rmplanet 
      write(5,26) sunpot,bsun  
  25  format(1x,'Radius Sun =',e11.3,' m, radius planet =',e11.3,
     1	' m, mass sun =',e11.3,' kg,   mass planet =', e11.3,' kg'/)
  26  format(1x,'Potential Sun =',e11.3,' V  Sun magnetic B field',
     1    e11.3,' T'/)
  27  format(1x,'    n        t           H        sx(m)       ',
     1    'sy(m)     velx(m/s) ',
     1'  vely(m/s)     gacc       eacc        maccx        maccy  '/)
  28  format(1x,i6,9e12.3)
  29  format(1x,'q Sun =',e11.3,' C  q Planet',e11.3,' C'/)
c-----------------------------------------------------------------
c g is gravitational constant; N m2 kg-2
c epsilon electrostatic constant; C V-1 m-1
      g = 6.67e-11
      pie = 3.14
      epsilon = 8.85e-12
      qsun = 4*pie*sunpot*epsilon*rsun
      qplanet = 4*pie*(2*sunpot/3)*epsilon*rplanet
      write(*,29) qsun,qplanet  
      write(5,29) qsun,qplanet
      H = dt*86400
c dt is in days, hence factor 86400 to bring to seconds
c H,HMAX,HMIN are step sizes, CREL,CABS, relative and absolue errors.
      HMAX = 10000*H
      HMIN = H/10000
      CREL = 0.001
c      CABS = 1.0E-10
      CABS = 1.0E9
      CFACT = 5  
      t = 0.0
      velx = velplanet
      vely = 0.0
      velz = 0.0
      Y(1) = sx
      Y(2) = velx
      Y(3) = sy
      Y(4) = vely
      Y(5) = sz
      Y(6) = velz	
      nprint = ndprint
      write(*,27)
      write(5,27)
      n = 0
      write(*,28) n,t,H,sx,sy,velx,vely
      write(5,28) n,t,H,sx,sy,velx,vely
      call HAMMING(6,t,H,HMAX,HMIN,CREL,CABS,CFACT,Y,F,OUTPUT,
     1    FAIL)
      close(5)
      pause
      end 
c  end main
c--------------------------------------------------------------------------
c      subroutine dipole(x,y,z,bx,by,bz)
c calculates magnetic fields bx,by,bz at any position x,y,z with the sun 
c at the origin assuming a dipole representation of the mnagnetic field such
c that the field at the surface of the sun at the sun's equatior is bsun.
c--------------------------------------------------------------------------
      subroutine dipole(x,y,z,bx,by,bz)
      common bsun,rsun,rplanet,g,epsilon,sunpot,pie,qsun,qplanet,
     1    rmplanet,rmsun,n,nsteps,ndprint,nprint
c dipole position at 0,0,zd and 0,0,-zd at 2/3 of sun's radius
      zd = 2*rsun/3
c      rmag = bsun*(rsun**2 + zd**2)**1.5/(2*zd)
      r1squared = x**2 + y**2 + (z-zd)**2
      r1 = sqrt(r1squared)
      r2squared = x**2 + y**2 + (z+zd)**2
      r2 = sqrt(r2squared)
      rxy = sqrt(x**2 + y**2)
      b1 = (bsun/r1squared)*(rsun**2 + zd**2)**1.5/(2*zd)
      b2 = -(bsun/r2squared)*(rsun**2 + zd**2)**1.5/(2*zd)
c      rmag = bsun*(rsun**2 + zd**2)**1.5/(2*zd)
c      b1 = rmag/r1squared
c      b2 = -rmag/r2squared
c b1x,b1y,b1z are components of first dipole force in x, y and z direction.
      b1x = b1*x/rxy
      b1y = b1*y/rxy
      b1z = b1*(z-zd)/r1
      b2x = b2*x/rxy
      b2y = b2*y/rxy
      b2z = b2*(z+zd)/r2
      bx = b1x + b2x
      by = b1y + b2y
      bz = b1z + b2z
  28  format(1x,i6,6e12.3)
      end
c  end subroutine dipole
C--------------------------------------------------------------------------
  
      SUBROUTINE F(XX,Y,DY,IC)
      common bsun,rsun,rplanet,g,epsilon,sunpot,pie,qsun,qplanet,
     1    rmplanet,rmsun,n,nsteps,ndprint,nprint
      common/out/ sx,sy,velx,vely,gacc,eacc,rmaccx,rmaccy,rmaccz,H,s
      DIMENSION Y(6),DY(6)
c s is distance from center of sun
c sxy is distance projected on xy plane
c accelleration due to gravity; gacc
      sx = Y(1)
      velx = Y(2)
      sy = Y(3)
      vely = Y(4)
      sz = Y(5)
      velz = Y(6)
      s = sqrt(sx**2 + sy**2 +sz**2)
      sxy = sqrt(sx**2 + sy**2) 
      gacc = -g*rmsun/s**2
      gaccx = gacc*sx/sxy
      gaccy = gacc*sy/sxy
      gaccz = gacc*sz/s
c accelleration due to electrostatic force; eacc
      eacc = qsun*qplanet/(4*pie*epsilon*rmplanet*s**2)
      eaccx = eacc*sx/sxy
      eaccy = eacc*sy/sxy
      eaccz = eacc*sz/s
c accelleration due to magnetic field; rmacc
      call dipole(sx,sy,sz,bx,by,bz)
      rjx = qplanet*velx
      rjy = qplanet*vely
      rjz = qplanet*velz
      rmaccx = (rjy*bz - rjz*by)/rmplanet
      rmaccy = -(rjx*bz - rjz*bx)/rmplanet
      rmaccz = (rjx*by - rjy*bx)/rmplanet
c add accelleration components
      accx = gaccx + eaccx + rmaccx
      accy = gaccy + eaccy + rmaccy
      accz = gaccz + eaccz + rmaccz
c Y(1) = sx, Y(2) = velx, Y(3) = sy, Y(4) = vely, Y(5) = sz, Y(6) = velz
c DY(1) = velx, DY(2) = -accx, DY(3) = vely, DY(4) = -accy, DY(5) = velz, DY(6) = -accz
      DY(1) = Y(2)
      DY(2) = accx
      DY(3) = Y(4)
      DY(4) = accy
      DY(5) = Y(6)
      DY(6) = accz
      RETURN
      END
c---------------------------------------------------------------------------------------
		  
      SUBROUTINE OUTPUT(XX,Y,DY,IC)
      common bsun,rsun,rplanet,g,epsilon,sunpot,pie,qsun,qplanet,
     1    rmplanet,rmsun,n,nsteps,ndprint,nprint
      common/out/ sx,sy,velx,vely,gacc,eacc,rmaccx,rmaccy,rmaccz,H,s
      DIMENSION Y(6),DY(6)
  28  format(1x,i6,10e12.3)
      n = n+1
      if (n.eq.nprint) then
	  write(*,28) n,XX,H,sx,sy,velx,vely,gacc,eacc,rmaccx,rmaccy
	  write(5,28) n,XX,H,sx,sy,velx,vely,gacc,eacc,rmaccx,rmaccy
	  nprint = nprint + ndprint
      endif
      if (s.lt.rsun) IC = 1
      if (n.ge.nsteps) IC = 1
  2   RETURN
      END    		  
c----------------------------------------------------------------

      SUBROUTINE FAIL(XX,Y,DY,H)
          DIMENSION Y(6),DY(6)
    1     FORMAT (F6.3,3E12.4)
    2     FORMAT ('   FAILED AT')
          WRITE (5,2)
          WRITE (5,1) XX,Y(1),Y(2),H
          WRITE (*,2)
          WRITE (*,1) XX,Y(1),Y(2),H
          RETURN
          END
C----------------------------------------------------------------

      SUBROUTINE HAMMING(K,X,H,HMAX,HMIN,CREL,CABS,ICFACT,YC,
     1    BOXA,BOXB,BOXD)
C HAMMING METHOD FOR SOLVING A SYSTEM OF FIRST ORDER ORDINARY
C     DIFFERENTIAL EQUATIONS
C     K  NUMBER OF EQUATIONS
C     X  INDEPENDENT VARIABLE CONTAINING ITS STARTING VALUE
C     H  VARIABLE CONTAINING STARTING VALUE OF STEP SIZE
C     HMAX MAXIMUM ALLOWABLE STEP SIZE
C     HMIN  MINIMUM ALLOWABLE STEP SIZE
C     CREL  RELATIVE ERROR TOLERANCE
C     CABS  ABSOLUTE ERROR TOLERANCE
C     CFACT  FACTOR USED FOR THE CRITERION FOR DOUBLING H
C     YO ARRAY OF INITIAL VALUES OF DEPENDENCY Y VARIABLES
C     BOXA PROCEDURE FOR EVALUATING THE DERIVATIVES 
C     BOXB PROCEDURE FOR OUTPUT AFTER EACH SUCCESSFUL INTEGRATION   
C     BOXD  PROCEDURE FOR EXITING WHEN ERROR CRITERIA FAILS
C
      INTEGER K,ICFACT
      DIMENSION Y(7,9),DY(7,9),C(7),DC(7),YY(7),DYY(7),RK1(7),
     1 RK2(7),RK3(7),DTABS(7),YC(7)
      INTEGER IA,IB,I,IM1,J,L,LP1,LP2,LP3,LP4,I2M1
      EXTERNAL BOXA,BOXB,BOXD
      IC = 0
      RTEST=15*CREL
      ATEST=15*CABS
      DO 10 J=1,K
  10      Y(J,1)=YC(J)
          IF(HMAX.LT.0) H=-ABS(H)
C  THE NEXT TWO STATEMENTS ADDED BY RWL FOR OBTAINING INITIAL CONDITIONS 2/5/75
          CALL BOXA(X,YC,DYY,IC)    
          IF (IC.EQ.1) RETURN
          CALL BOXB(X,YC,DYY,IC)
          IF (IC.EQ.1) RETURN
C RUNGE-KUTTA FOURTH ORDER FOR STARTING INTEGRATION
          H=2*H
          IA=2
          IB=2
 100  DO 190 I=IA,IB
          IM1=I-1
      DO 110 J=1,K
  110     YY(J)=Y(J,IM1)
          CALL BOXA(X,YY,DYY,IC)
          IF (IC.EQ.1) RETURN
      DO 120 J=1,K
          T=DYY(J)
          DY(J,IM1)=T
          T=T*H
          RK1(J)=T
          T=YY(J)+0.5*T
          YY(J)=T
          Y(J,I)=T
 120  CONTINUE
          XTEMP=X+.5*H
          CALL BOXA(XTEMP,YY,DYY,IC)
          IF (IC.EQ.1) RETURN
      DO 130 J=1,K
         T=DYY(J)
         DY(J,I)=T
         T=T*H
         RK2(J)=T
         T=Y(J,IM1)+.5*T
         YY(J)=T
         Y(J,I)=T
 130  CONTINUE
      CALL BOXA(XTEMP,YY,DYY,IC)
      IF (IC.EQ.1) RETURN
      DO 140 J=1,K
          T=DYY(J)
          DY(J,I)=T
          T=T*H 
          RK3(J)=T
          T=Y(J,IM1)+T
          YY(J)=T
          Y(J,I)=T
 140  CONTINUE
          X=X+H	
          CALL BOXA(X,YY,DYY,IC)
          IF (IC.EQ.1) RETURN
      DO 150 J=1,K
         Y(J,I)=Y(J,IM1)+.166666666666*(RK1(J)+2*(RK2(J)+RK3(J))+
     1    H*DYY(J))
  150 	  DY(J,I)=DYY(J)
  190 CONTINUE
  200 IF(.NOT.(IB.EQ.2)) GO TO 250
  210 DO 220 J=1,K
  220     C(J)=Y(J,2)
      X=X-H
      H=H/2.
      IF(ABS(H).GE.ABS(HMIN)) GOTO 240
      DO 230 J=1,K
          YY(J)=Y(J,1)
 230  DYY(J)=DY(J,I)
      CALL BOXD(X,YY,DYY,H)
      RETURN
 240  IB=3
      GOTO 100
 250  IF(.NOT.(IB.EQ.3)) GOTO 300
      DO 260 J=1,K
          T=Y(J,3)
          TEST=ABS(C(J)-T)
          TT=ABS(T)*RTEST
          IF(.NOT.((TEST.GT.TT).AND.(TEST.GT.ATEST))) GOTO 260
          X=X-H
          GOTO 210
 260  CONTINUE
      IA=4
      IB=4
          GO TO 100
 300  CONTINUE
      DO 310 J=1,K
 310      YY(J)=Y(J,4)
          CALL BOXA(X,YY,DYY,IC)
          IF (IC.EQ.1) RETURN
      DO 320 J=1,K
 320      DY(J,4)=DYY(J)
C OUTPUT RUNGE-KUTTA RESULTS
          X=X-4.*H
      DO 400 I=1,4
          X=X+H
      DO 330 J=1,K
	  YY(J)=Y(J,I)
 330      DYY(J)=DY(J,I)
 400  CONTINUE
      CALL BOXB(X,YY,DYY,IC)
      IF (IC.EQ.1) RETURN
C CONSTANTS AND TEST VARIABLES FOR HAMMING
      CON2=1./121.
      CON1=112.*CON2
      CON2=9.*CON2
      CON3=1./128.
      T=1./CON2
      TREL=CREL*T
      TABS=CABS*T
      DTABS(1)=TABS*1.0E-4
      DTABS(2)=TABS*1.0E-5
      DTABS(3)=TABS*1.0E-2
      T=ICFACT
      IF(T.EQ.0.0) T=100.
      TDH=TREL/T
      L=1
      GOTO 510
C RESTORE Y AND DY
 410  DO 420 I=1,4
          LPI=L+I
      DO 420 J=1,K
          Y(J,I)=Y(J,LPI)
          DY(J,I)=DY(J,LPI)
 420  CONTINUE
      L=1
C INTEGRATE
 510  X=X+H
      LP1=L+1
      LP2=L+2
      LP3=L+3
      LP4=L+4
      DO 520 I=1,K
          C(I)=Y(I,L)+1.333333*H*(2*(DY(I,LP1)+DY(I,LP3))-
     1    DY(I,LP2))
 520  CONTINUE
      CALL BOXA(X,C,DC,IC)
      IF (IC.EQ.1) RETURN
      DO 600 I=1,K
          T=.125*(9*Y(I,LP3)-Y(I,LP1)+3*H*(DC(I)+2*DY(I,LP3)-
     1    DY(I,LP2)))
          C(I)=C(I)-T
          Y(I,LP4)=T
          TT=TREL*ABS(T)
          TEST=ABS(C(I))
          TABS=DTABS(I)
          IF((TEST.GT.TT).AND.(TEST.GT.TABS)) GOTO 710
 600  CONTINUE
 610  DO 620 I=1,K
          YY(I)=Y(I,LP4)+CON2*C(I)
 620      Y(I,LP4)=YY(I)
      CALL BOXA(X,YY,DYY,IC)
      IF (IC.EQ.1) RETURN
      CALL BOXB(X,YY,DYY,IC)
      IF(IC.EQ.1) RETURN
      DO 630 I=1,K
 630      DY(I,LP4)=DYY(I)
C CHECK FOR DOUBLING INTEGRATION STEP
          IF(ABS(2.*H).GT.ABS(HMAX)) GOTO 410
      DO 640 I=1,K
          TT=TDH*ABS(YY(I))
          TEST=ABS(C(I))
          TABS=DTABS(I)
          IF((TEST.GT.TT).AND.(TEST.GT.TABS)) GOTO 410
 640  CONTINUE
	  IF(.NOT.L.EQ.3) GOTO 700
      DO 650 I=1,4
          I2M1=2*I-1
      DO 650 J=1,K
          Y(J,I)=Y(J,I2M1)
          DY(J,I)=DY(J,I2M1)
 650  CONTINUE
      H=2*H
      L=1
      GO TO 510
 700  L=LP1
      GO TO 510
C CHECK FOR HALVING STEP SIZE. IF LESS THAN HMIN, TERMINATE
 710  IF(ABS(.5*H).GE.ABS(HMIN)) GOTO 730
      DO 720 I=1,K
          YY(I)=Y(I,LP3)
 720      DYY(I)=DY(I,LP3)
      X=X-H
      CALL BOXD(X,YY,DYY,H)
      RETURN
 730  DO 740 J=1,K
          Y1=Y(J,LP1)
          Y2=Y(J,LP2)
          Y3=Y(J,LP3)
          DY1=DY(J,LP1)
          DY2=DY(J,LP2)
          DY3=DY(J,LP3)
          Y(J,1)=CON3*(11*Y3+72*Y2+45*Y1-H*(3*DY3+36*DY2-
     1        9*DY1))
          Y(J,2)=Y2
	  DY(J,2)=DY2
          Y(J,3)=CON3*(45*Y3+72*Y2+11*Y1+H*(-9*DY3+36*DY2+
     1        3.*DY1))
          Y(J,4)=Y3
          DY(J,4)=DY3
 740  CONTINUE
      X=X-H
      H=H/2.
      L=1
      T=X-3.*H
      TT=X-H
      DO 750 I=1,K
          C(I)=Y(I,1)
 750      YY(I)=Y(I,3)
      CALL BOXA(T,C,DC,IC)
      IF (IC.EQ.1) RETURN
      CALL BOXA(TT,YY,DYY,IC)
      IF (IC.EQ.1) RETURN
      DO 760 I=1,K
          DY(I,1)=DC(I)
 760      DY(I,3)=DYY(I)
      GO TO 510
C END HAMMING
      END
