C C PURPOSE: CALCULATES A MAGNETIC FIELD ACCORDING TO ZHANG AND LOW C C INPUTS: C RAD RADIUS OF POINT RELATIVE TO SUN CENTER, UNITS OF SOLAR RADII C THET AZIMUTHAL ANGLE, UNITS RADIANS, IN S FRAME C PH POLAR ANGLE, UNITS RADIANS, IN S FRAME C IPRINT SET TO 1 FOR MORE VERBOSE OUTPUT C OUTPUTS: C B MAGNITUDE OF FIELD C T POLAR ANGLE C P AZIMUTHAL ANGLE C EDENS ELECTRON DENSITY /CM3 C HDENS HYDROGEN NUCLEI DENSITY /CM3 C TEMP ELECTRON TEMPERATURE K C VEL FLOW VELOCITY ALONG LINE OF SIGHT (+VE MEANS TOWARDS OBSERVER) C TURBV MICROTURBULENCE KM/S (CONTRIBUTES TO LINE PROFILE WIDTHS) C COMMON: C C COMMENTS: APRIL 19, 2006 C C NOTES: USER-SUPPLIED ROUTINE C SUBROUTINE ZL(RAD,THET,PH,BMAG,T,P,EDENS,HDENS, * TEMP,VEL,TURBV,IPRINT) INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CCONST' INCLUDE 'CGRID' INCLUDE 'CINPUT' INCLUDE 'CLU' INCLUDE 'CINTS' COMMON /CZL/ ACONST,B0,CI0, ALPHAX, ALPHAY,P0,T0 COMMON /WATM/BX,BY,BZ DATA KCALL /0/ DATA EPSILON, ZER, SMALL /1.E-3,0.,1.E-20/ SAVE KCALL C IF (KCALL .EQ. 0) THEN CALL CSTRIP(LDIP,'ZL.DAT') READ(LDUMS,*) ACONST READ(LDUMS,*) B0 READ(LDUMS,*) CI0 READ(LDUMS,*) ALPHAY, ALPHAX READ(LDUMS,*) P0 READ(LDUMS,*) T0 WRITE(LOUT,4001) ACONST,B0,CI0,ALPHAY,ALPHAX,P0,T0 4001 FORMAT('1 CORONAL ROUTINE IS ZL'/ * 'A= ',F11.6,' B0= ',F10.4,' CI0= ',F10.4,' ALPHAY= ',F10.4, * ' ALPHAX= ',F10.4/'P0= ',F10.4,' T0=',1P,E10.3,0P) CALL CLOSE(LDUMS) KCALL=KCALL+1 ENDIF C C CALCULATE CARTESIAN COORDINATES OF THE POINT RAD,THET,PH C X IS LOS C X=RAD*SIN(THET)*COS(PH) Y=RAD*SIN(THET)*SIN(PH) Z=RAD*COS(THET) R=RAD C C ROTATE TO (R,THETA) PLANE IN WHICH FFL SOLUTION IS GIVEN IN TERMS C OF STREAM FUNCTION. C C ROTATION ABOUT Y AXIS, AND THEN X-AXIS C X1=X Y1=Y Z1=Z CALL ROT(X1,Z1,ALPHAY) CALL ROT(Y1,Z1,ALPHAX) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ROTATED FRAME: C CALCULATE THE (R,THETA) COMPONENTS OF THE MAGNETIC FIELD IN THE ROTATED C REFERENCE FRAME. C COSTH=Z1/R SINTH=SQRT(1.-COSTH*COSTH) TH=ATAN2(SINTH,COSTH) TANTH=SINTH/COSTH C DTH=0.01 DR=0.001*R C C SIMPLE CENTERED DIFFERENCE, 4 CALLS C C BR=1./R/R/SINTH*(AZL(R,TH+DTH,CI0,ACONST,B0) - C * AZL(R,TH-DTH,CI0,ACONST,B0))/2./DTH C BTH=-1./R/SINTH * (AZL(R+DR,TH,CI0,ACONST,B0) - C * AZL(R-DR,TH,CI0,ACONST,B0))/2./DR C ALL DONE! C C SIMPLE CENTERED DIFFERENCE, NR ITERATION C BR=1./R/R/SINTH*DAZLDTH(R,TH,CI0,ACONST,B0) BTH=-1./R/SINTH*DAZLDR(R,TH,CI0,ACONST,B0) C C CARTESIAN COMPONENTS IN THIS ROTATED REFERENCE FRAME C BMAG=SQRT(BR*BR+BTH*BTH) V=SQRT(X1*X1+Y1*Y1) COSPH=X1/V SINPH=Y1/V BZ= BR*COSTH - BTH*SINTH BX= (BR*SINTH + BTH*COSTH)*COSPH BY= (BR*SINTH + BTH*COSTH)*SINPH C C DENSITY RESTRICTION AROUND CURRENT SHEET, NARROWING WITH R C NOT IN PRESSURE EQUILIBRIUM C BEGIN C THICK=0.1/R/R C TH1=ACOS(COSTH) C DELTA=(PI/2.-TH1)/THICK C END C ROTATED FRAME END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TRANSFORM BACK, ROTATE ABOUT X AND Y AXES C BX1=BX BY1=BY BZ1=BZ CALL ROT(BY1,BZ1,-ALPHAX) CALL ROT(BX1,BZ1,-ALPHAY) C RE-ASSIGN BX=BX1 BY=BY1 BZ=BZ1 C C RETURN VALUES OF B,T,P C BMAG=SQRT(BX*BX+BY*BY+BZ*BZ) T=ACOS(BZ/BMAG) P=ATAN2(BY,BX) C C THERMAL PARAMETERS C C C ELECTRON DENSITY AND TEMPERATURE. C VEL=0. TURBV=15. C DELTA=0. C GFAC=EXP(-DELTA*DELTA) GRAV=2.7398E4 TEMP=T0 C EDENS=P0/BK/T0/1.88*GFAC EDENS=P0/BK/T0/1.88 RHO=1.3625*0.8*UU PTOT=BK*TEMP*1.88 + 0.5* RHO*(TURBV*TURBV)*1.E10 HSCALE = PTOT/RHO/GRAV HSCALE=HSCALE/RSUNCM FFAC= MAX((R-1.),ZER) FFAC=EXP(-FFAC/HSCALE) EDENS=EDENS*FFAC HDENS=0.8*EDENS RETURN END C C********************************************************************** C FUNCTION AZL(R, THETA, CI0, A0, B0) INCLUDE 'PREC' INCLUDE 'CCONST' EPS = 1.E-6 AZL=0.0 ST=SIN(THETA) C C ABRAMOWITZ & STEGUN EQS. 25.4.38, INTEGRAL OF F(X)/SQRT(1-X^2) DX BETWEEN C -1 AND 1 C HERE, F(X) IS X/SQRT(A^2+R^2 - 2AR SIN(TH) SQRT(1-X^2)) C NMAX=40 SUM=0. W = PI/NMAX DO 100 I=1,NMAX X = COS( PI*(2.*I-1.)/NMAX ) F = X/SQRT(EPS+A0*A0+R*R -2.*A0*R*ST*X) FIMAGE = X/SQRT(EPS+A0*A0*R*R +1. -2.*A0*R*ST*X) SUM=SUM+F-FIMAGE 100 CONTINUE AZL = 2.*CI0*A0*R*ST*SUM*W + B0*ST*ST/R RETURN END C C********************************************************************** C BASED ON DFRIDR.F OF NUMERICAL RECIPES C FUNCTION DAZLDTH(R,TH,CI0,ACONST,B0) INTEGER NTAB REAL DAZLDTH,ERR,H,X,FUNC,CON,CON2,BIG,SAFE PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.) INTEGER I,J REAL ERRT,FAC,HH,A(NTAB,NTAB) X=TH H=0.01 HH=H A(1,1)=(AZL(R,TH+HH,CI0,ACONST,B0)-AZL(R,TH-HH,CI0,ACONST,B0)) * /(2.0*HH) ERR=BIG DO 12 I=2,NTAB HH=HH/CON A(1,I)=(AZL(R,TH+HH,CI0,ACONST,B0)-AZL(R,TH-HH,CI0,ACONST,B0)) * /(2.0*HH) FAC=CON2 DO 11 J=2,I A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.) FAC=CON2*FAC ERRT=MAX(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1))) IF (ERRT.LE.ERR) THEN ERR=ERRT DAZLDTH=A(J,I) ENDIF 11 CONTINUE IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR) RETURN 12 CONTINUE RETURN END C C********************************************************************** C FUNCTION DAZLDR(R,TH,CI0,ACONST,B0) INTEGER NTAB REAL DAZLDTH,ERR,H,X,FUNC,CON,CON2,BIG,SAFE PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.) INTEGER I,J REAL ERRT,FAC,HH,A(NTAB,NTAB) X=R H=0.01*R HH=H A(1,1)=(AZL(R+HH,TH,CI0,ACONST,B0)-AZL(R-HH,TH,CI0,ACONST,B0)) * /(2.0*HH) ERR=BIG DO 12 I=2,NTAB HH=HH/CON A(1,I)=(AZL(R+HH,TH,CI0,ACONST,B0)-AZL(R-HH,TH,CI0,ACONST,B0)) * /(2.0*HH) FAC=CON2 DO 11 J=2,I A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.) FAC=CON2*FAC ERRT=MAX(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1))) IF (ERRT.LE.ERR) THEN ERR=ERRT DAZLDR=A(J,I) ENDIF 11 CONTINUE IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR)RETURN 12 CONTINUE RETURN END