module coef_module ! use param_module,only: rmlat,rmlon ! implicit none ! integer,parameter :: | nm = 80, | mm = 24, | nlm = rmlon, | nthp = rmlat + 1, | tmpo = 2*mm + 1, | nmp = nm + 1, | mmp = mm + 1, | ny = nlm, | nth = rmlat, | nqmax = 9, ! for qnm (qost in code) max 9 if less comment out qost | mqmax = 9, ! for qnm (qost in code) max 9 if less comment out qost | iodd_qnm = 1, ! take odd (n-m) for q_nm into accounts if 0 comment out qost | icm3e = 0 ! relate internal & external comp. by cm3e NASA TM 2000 ! eq(15-16) real :: qost(0:nqmax,0:mqmax), | c1(0:nm,-mm:mm), | c2(0:nqmax,-mqmax:mqmax), | r(0:nm,0:mm), ! common/eqcr/ | pmopmmo(0:mm), ! common/eqcr/ | p(0:nm,0:mm), ! spherical harmonics | dp(0:nm,0:mm), ! derivative of spherical harmonics | f(-mm:mm), ! common/far/ | c_tor(0:nm,-mm:mm), ! for calculating the toroidal eqc | c_fld(0:nm,-mm:mm), ! for calculating the fld-algnd eqc | c_tot(0:nm,-mm:mm), ! for calculating the total eqc | c_grd(0:nm,-mm:mm), ! for calculating EQC @ ground | c_mptext(0:nm,-mm:mm), ! magnetic potential external part | c_mptextgrd(0:nm,-mm:mm), ! magnetic potential external part @ ground | c_mptintgrd(0:nm,-mm:mm), ! magnetic potential induced part @ ground | c_mpttotgrd(0:nm,-mm:mm) ! magnetic potential total (external+induced) @ ground complex :: | qnnmm(0:nm,0:mm) ! correlates external to internal coeff. from cm3e ! contains !----------------------------------------------------------------------- ! calulated the coefficients for the equivalent current by ! using the kqlam and kqphi ! couldn't represent the electrojet very good ! subroutine coefficient(h_lay) use param_module,only: re,qcnst1,qcnst2,sq2,rad,pi,rmlat use eqvcur_module,only: reglat,reglon,kqlam,kqphi implicit none real,intent(in) :: h_lay ! height of reference layer [km] real :: xmlon, date, colat,elon common/ssolar/ xmlon,date,colat,elon integer :: mp,nlmm,jmx,jcon,np,mq,mqf,mc,ias, | j,l,n,m,ii,jj,i, | in,im real :: xms,xns,den,rdphofp,dph,ctm,th,ctp,daofpr,ct,st,sts, | pm2,xl,x,s,ph,r_lay, | jsc,jec,jss,jsa,jes,jea,jph,jth ! r_lay = re + 1000.*h_lay ! radius reference layer [m] ! do MP = 1, MMP ! m+1 = 1:25 M = MP - 1 ! 0:24 XMS = float(M*M) ! m^2 IF (MP.NE.1) PMOPMMO(M) = SQRT(1. + .5/FLOAT(M)) ! sqrt(1+1/(2m)) do N=M,NM ! n = m,49 XNS = float(n*n) ! n^2 DEN = AMAX1(4.*XNS - 1.,1.) ! 4n^2-1 R(N,M) = SQRT((XNS - XMS)/DEN) ! R_n^m = sqrt[(n^2-m^2)/(4n^2-1)] C1(N,-M) = 0. ! initialize C1(N,+M) = 0. enddo enddo ! c2(:,:) = 0. ! initialize c_fld(:,:) = 0. qost(:,:) = 0. c_tor = 0. c_fld = 0. c_tot = 0. c_grd = 0. c_mptext = 0. c_mptextgrd = 0. c_mptintgrd = 0. c_mpttotgrd = 0. ! NLMM = NLM - 1 ! rmlon -1 JMX = NTHP/2 ! (rmlat+1)/2 RDPHOFP = .5*R_lay/NLMM !r_i/2/(rmlon-1) DPH = 2.*pi/NLMM ! 2Pi/(rmlon-1) CTM=1. ! cos(th) reglat(rmlat) = 88.98 ! to avoid numerical problems with +/- 90. reglat(1) = -88.98 ! latitude loop: J over N hemispheresince jcon is used DO 200 J=1,JMX ! 1:(rmlat+1)/2 JCON = rmlat - J + 1 ! index of conjugate point rmlat:(rmlat+1)/2 ! TH is colatitude (N pole = 0) th = 90. - 0.5*(reglat(jcon)+reglat(jcon-1)) th = th/rad ! deg to rad ctp = cos(th) ! cos(colatitude) IF (J.EQ.JMX) CTP = 0. ! DAOFPR is surface area element of integration, divided by 4pi*RI DAOFPR = RDPHOFP*(CTM - CTP) ! r_i/2/(rmlon-1)*delta(cos(th)) CTM = CTP th = 90. - reglat(jcon) ! colatitude th = th/rad ! deg to rad CT = COS(TH) ! cos(th) ST = SIN(TH) ! sin(th) STS = ST*ST ! sin^2(th) ! In this subroutine P = P_n^m/sin(th), where P_n^m is an associated Legendre ! polynomial, fully normalized as in Richmond (1974). ! e.g., P(0,0) = P(n=0,m=0)/sin(th). ! DP(n,m) = dP_n^m/d th ! P(0,0) = 1./st ! P(n=0,m=0)/sin(th) P(1,1) = PMOPMMO(1) ! P(n=0,m=0)/sin(th)*sin(th)*sqrt(1+1/(2m)) DO 80 MP=1,MMP ! m+1 = 1:25 M = MP - 1 ! 0:24 IF (M.GE.2) P(M,M) = PMOPMMO(M)*P(M-1,M-1)*ST ! P_n-1^m-1/sin(th)*sin(th)*sqrt(1+1/(2m)) DP(M,M) = M*CT*P(M,M) ! dP_n^m/d th = n*cos(th)*P_n^m/sin(th) PM2 = 0. DO 80 N=MP,NM ! m+1:nm (48) NP = N + 1 ! n+1 P(N,M) = (CT*P(N-1,M) - R(N-1,M)*PM2)/R(N,M) ! P_n^m/sin(th) = [cos(th)*P_n-1^m/sin(th) - ! R_n-1^m*P_n-2^m/sin(th)]/R_n^m PM2 = P(N-1,M) ! P_n-1^m/sin(th) DP(N,M) = N*CT*P(N,M) - (2*N + 1)*R(N,M)*PM2 ! dP_n^m/d th = n*cos(th)*P_n^m/sin(th)/cos(th) ! -(2n-1)*R_n^m*P_n-1^m/sin(th) 80 continue ! calculated Q_n^m/sin(th) paper Richmond [1974] (table 1 and eq. 30 from MATLAB code) ! parameter iodd_qnm = 0 take odd (n-m) for q_nm into accounts ! antisymmetric part: n-m even ! symmetric part: n-m odd ! call calqnm(ct,st,sts,th) DO 150 L=1,NLMM ! 1:rmlon-1 change from mag.long. to mag.local time PH = (reglon(l)+180.-xmlon)/rad ! mag.longitude - subsolar point [deg to rad] CALL FCMP(PH,MM) ! f_-m = ie^-(i m ph); f_m = ie^(i m ph) JSC = -kqlam(L,JCON) ! JSC southward current, N hemisphere JEC = kqphi(L,JCON) ! eastward current, N hemisphere JSS = JSC + kqlam(L,J) ! 2x symmetric component of equatorward current JSA = JSC - kqlam(L,J) ! 2x antisymmetric component of equatorward current JES = JEC + kqphi(L,J) ! 2x symmetric component of eastward current JEA = JEC - kqphi(L,J) ! 2x antisymmetric component of eastward current DO 120 m=-mm,mm ! loop over order m MP = 1 + IABS(M) ! abs(m) + 1 ! MC = MMP - M IAS = -1 ! asymmetric DO 120 NP=MP,NMP ! loop over degree n abs(m)+1:nm+1 N = NP - 1 ! n IAS = - IAS ! 1 -> abs(m)+1, abs(m)+3 ... ! IAS = 1 for spherical harmonics symmetric about the equator (related ! to equivalent current function for ANTISYMMETRIC currents); ! = -1 for spherical harmonics antisymmetric about the equator (related ! to equivalent current function for SYMMETRIC currents); IF (IAS.EQ.1) GO TO 105 !spherical harmonics symmetric JPH = JES ! symmetric component of eqc JTH = JSS GO TO 110 105 JPH = JEA ! antisymmetric of eqc JTH = JSA IF (NP.GT.nqmax+1) GO TO 110 ! no c_2 calculated ! ! Perform integrals over Earth for each spherical harmonic component ! Note: Table 1 of Richmond (1974) has a sign error for negative m: ! when m is negative, functions Qnm should be multiplied by -1. S = 1. IF (M.LT.0) S= -1. C2(N,M) = C2(N,M) + DAOFPR*((M*P(N,MP-1) !d z_n^m / d th -> [m P_n^m-Q_n^m/sin(th)*2*n*cos(th)]*f_-m*K_lam | - QOST(N,MP-1)*2.*CT*S*float(N))*F(-m)*JTH + | S*float(M)*QOST(N,MP-1)*F(M)*JPH) !d z_n^m / d ph -> m*f_m*Q_n^m/sin(th)*K_ph 110 C1(N,M) = C1(N,M)+ DAOFPR*(DP(N,MP-1)*F(M)*JPH ! grad Y_n^m = d P_n^m /d th f_m * K_ph - m P_n^m f_-m * K_th | - float(M)*P(N,MP-1)*F(-m)*JTH) 120 CONTINUE 150 CONTINUE 200 CONTINUE DO 250 m = -mm,mm MP = iabs(M)+1 DO 250 NP=MP,NMP N = NP - 1 den = N*(N+1) IF (DEN.EQ.0.) DEN = 1. C1(N,M ) = C1(N,M)/DEN ! * 1/ [n(n+1)] IF (N.LE.nqmax) C2(N,M) = C2(N,M)/DEN ! * 1/ [n(n+1)] 250 CONTINUE c_tor = c1 ! toroidal component c_tot = c1 ! total component ! calculate the coefficients for the total equivalent current function DO 300 m = -mqmax,mqmax MP = IABS(m) + 1 DO 300 N=MP, nqmax+1 c_tot(n-1,m) = c1(n-1,m) + c2(n-1,m) ! total c_fld(n-1,m) = c2(n-1,m) ! field-aligned 300 continue return end subroutine coefficient !---------------------------------------------------------------------------- ! calculates F(+ M) = SQ2*sin(m*ph) ! F(- M) = SQ2*cos(m*ph) ! SUBROUTINE FCMP(PH,MM) use param_module,only: sq2 implicit none integer,intent(in) :: mm real,intent(in) :: ph integer :: m,mmo,jj real :: cp,sp F(0) = 1. CP = COS(PH) SP = SIN(PH) F(- 1) = SQ2*CP F(+ 1) = SQ2*SP DO 85 M=2,MM MMO = M - 1 F(- M) = CP*F(- MMO) - SP*F(+ MMO) F(+ M) = CP*F(+ MMO) + SP*F(- MMO) 85 continue RETURN end subroutine fcmp !---------------------------------------------------------------------------- subroutine coef_eqcgrd ! original coefficients for calculating the influence from external field ! mapped from the ionospheric layer (110km) to the ground use param_module,only: re,ri,nue implicit none integer :: np,n,m,mp real :: aob,aobn,fak(0:nm) ! AOB = re/ri ! r(earth)/r(ionosphr) AOBN = ri/re ! do NP=1,NMP AOBN = AOBN*AOB N = NP - 1 FAK(N) = AOBN ! (re/ri)^(n) enddo ! c_grd = c_tot ! do m = -mm,mm ! c_grd = (re/ri)^(n) c_tot do n= iabs(m),nm c_grd(N,M ) = c_tot(N,M)*FAK(N) enddo enddo return end subroutine coef_eqcgrd !---------------------------------------------------------------------------- subroutine coef_magpot_ex ! calculate coefficient for mag. potential ! use param_module,only: nue ! implicit none ! integer :: np,n,m real :: fak(0:nm) ! ! v_nm coefficient for external magnetic potential (due to the ionosphere) DO NP=1,NMP ! nue (n+1)/(2n+1) N = NP - 1 FAK(N) = nue*float(NP)/FLOAT(2*N + 1) enddo ! c_mptext = c_grd ! c_mptext = nue (n+1)/(2n+1) c_grd do m = -mm,mm do n= iabs(m),nm c_mptext(N,M) = c_grd(N,M)*FAK(N) enddo enddo ! return end subroutine coef_magpot_ex !---------------------------------------------------------------------------- subroutine coef_magpot_tot ! calculate coefficient for total mag. potential (external + induced) ! use param_module,only: re,ri,depth ! re[m],ri[m],depth[km] implicit none ! ! v_nm coefficient for magnetic potential integer :: np,n,m real :: coa,coatnpo,coas,fak(nmp),cint COA = (re - DEPTH*1000.)/re ! mapping from re to re-depth COATNPO = 1./COA COAS = COA*COA do NP=1,NMP COATNPO = COATNPO*COAS N = NP - 1 FAK(N) = 1. + float(N)*COATNPO/FLOAT(NP) !1+ n/(n+1)[(r_e-depth)/r_e]^(2n+1) enddo if(icm3e == 0 ) then ! w/o cm3e relation internal/external c_mpttotgrd = c_mptext do m = -mm,mm do n= iabs(m),nm np = n + 1 c_mptextgrd(n,m) = c_mptext(N,M)/re ! c_mptext/re c_mptintgrd(n,m) = c_mptextgrd(N,M)*(FAK(N)-1.) ! c_mptext/re n/(n+1)[(r_e-depth)/r_e]^(2n+1) C_mpttotgrd(N,M) = c_mptext(N,M)*fak(n) ! c_mptext [1+ n/(n+1)[(r_e-depth)/r_e]^(2n+1)] enddo enddo elseif(icm3e == 1) then ! with cm3e relation internal/external c_mpttotgrd = c_mptext do m = -mm,mm do n= iabs(m),nm np = n + 1 c_mptextgrd(n,m) = c_mptext(N,M)/re ! c_mptext/re if(m /= 0) then ! c_mptintgrd(n,m) = c_mptextgrd(N,M)* ! re(q_n^m)*c_ext_n^m + im(q_n^m)*c_ext_n^-m ! | real(qnnmm(n,iabs(m))) +m/abs(m)*c_mptextgrd(N,-M)* ! | imag(qnnmm(n,iabs(m))) c_mptintgrd(n,m) = c_mptextgrd(N,M)* ! re(q_n^m)*c_ext_n^m + im(q_n^m)*c_ext_n^-m | real(qnnmm(n,iabs(m)))-m/abs(m)*c_mptextgrd(N,-M)* | imag(qnnmm(n,iabs(m))) ! write (6,'(1x,2i5,2f8.4)') N,m, ! * real(qnnmm(n,iabs(m))),imag(qnnmm(n,iabs(m))) else ! m = 0 with depth = 1000. km fak(n) = float(n)/float(n+1)*((re-1000.)/re)**(2*n+1) ! c_mptext/re n/(n+1)[(r_e-depth)/r_e]^(2n+1) c_mptintgrd(n,m) = c_mptextgrd(N,M)*fak(n) endif C_mpttotgrd(N,M) = re*(c_mptextgrd(N,M)+c_mptintgrd(n,m)) enddo enddo endif ! return end subroutine coef_magpot_tot !----------------------------------------------------------------------- subroutine calqnm(ct,st,sts,phi) ! calculate QOST = Qnm of Richmond (1974) eq.30 divided by sin(th) ! Index labeling is ! n , m , e.g., QOST(0,0) = Q(n=0,m=0)/sin(th). ! am 10/06/02 ! original goes upto 6 extended to higher order terms ! higher order terms are calculated with Mathematica4 on meeker ! $fshome/Mathematica/intlegp.nb use param_module,only: qcnst1,qcnst2 implicit none integer :: i real,intent(in) :: ct,st,sts,phi real :: x,xl,csc,cot qost(:,:) = 0. ! XL = 1. + .5*STS*ALOG((1. + CT)/(1. - CT))/CT ! L csc = 1./sin(phi) cot = cos(phi)/sin(phi) ! QOST(1,1) = 1.2247449*CT ! n=1 m=1 QOST(2,2) = 1.3693064*CT*ST*XL ! n=1 m=2 QOST(3,1) = .2291288*CT*(4. - STS*(3. + 6.*STS)) ! n=3 m=1 QOST(3,3) = 1.4790199*CT*STS*(1. + 2.*STS) ! n=3 m=3 X = STS*(2. + 3.*STS*XL) QOST(4,2) = QCNST1*CT*ST*(4. - X) ! n=4 m=2 QOST(4,4) = QCNST2*CT*ST*X ! n=4 m=4 X = STS*(3. + STS*(4. + 8.*STS)) QOST(5,1) = .018021728*CT*(56. - STS*(188. - 13.*X)) ! n=5 m=1 QOST(5,3) = .5255737*CT*STS*(8. - X) ! n=5 m=3 QOST(5,5) = .5484353*CT*STS*X ! n=5 m=5 X = STS*(8. + STS*(10. + 15.*STS*XL)) QOST(6,2) = .057727979*CT*ST*(64. - STS*(168. - 3.*X)) ! n=6 m=2 QOST(6,4) = .079047290*CT*ST*STS*(80. - 3.*X) ! n=6 m=4 QOST(6,6) = .2140611*CT*ST*STS*X ! n=6 m=6 qost(7,1) = 2.727689212715145*Cos(phi)*Sqrt(Sin(phi)**2)* ! n=7 m=1 - (-0.8807848156959635 + Sin(phi)**2)* - (-0.5697970475304803 + Sin(phi)**2)* - (-0.20824163331106132 + Sin(phi)**2)* - (1.6722005707476288 + Sin(phi)**2)* - (2.3384358827460394 + 0.4866229257898761*Sin(phi)**2 + - Sin(phi)**4) qost(7,3) = -6.5031302062421315*Cos(phi)*(Sin(phi)**2)**1.5* ! n=7 m=3 - (-0.8667204829110627 + Sin(phi)**2)* - (-0.5130377430396306 + Sin(phi)**2)* - (1.5075571030560535 + Sin(phi)**2)* - (1.9628371080697873 + 0.3722011228946397*Sin(phi)**2 + - Sin(phi)**4) qost(7,5) = 6.937218462755805*Cos(phi)*Sin(phi)**4* ! n=7 m=5 - Sqrt(Sin(phi)**2)*(-0.8233156743959129 + Sin(phi)**2)* - (1.1758902498268906 + Sin(phi)**2)* - (1.2911504175766504 + 0.14742542456902236*Sin(phi)**2 + - Sin(phi)**4) qost(7,7) =-0.0006924195703560433*Sqrt(Sin(phi)**2)* ! n=7 m=7 - (800.0000000000014*Cos(phi) - 1650.0000000000007*Cos(3*phi) + - 1289.9999999999995*Cos(5*phi) - 572.*Cos(7*phi) + - 156.00000000000003*Cos(9*phi) - 26.*Cos(11*phi) + - 2.*Cos(13*phi) + (0,2.2737367544323206e-13)*Sin(phi) + - (0,2.2737367544323206e-13)*Sin(3*phi) + (0,0.)*Sin(5*phi) + - (0,0.)*Sin(7*phi) + (0,0.)*Sin(9*phi) + - (0,0.)*Sin(11*phi) + (0,0.)*Sin(13*phi)) qost(7,:) = -1.*qost(7,:)/st qost(8,2) = Sin(phi)**2*(2.110733549333601*Cos(phi) + ! n=8 m=2 - 1.4532447593961175*Cos(3*phi) + - 1.4111626853260841*Cos(5*phi) + - 0.32273059308222324*Cos(7*phi) - - 0.08303591121779125*Cos(9*phi) + - 0.013114300307894611*Cos(11*phi) - - 0.0009595829493581424*Cos(13*phi) + - (-3.9304517605709512*Log(Cos(phi/2.)) + - 3.9304517605709512*Log(Sin(phi/2.)))*Sin(phi)**14) qost(8,4) = Sin(phi)**4*(8.932170688397838*Cos(phi) + ! n=8 m=4 - 6.2191283263968415*Cos(3*phi) + - 1.1286386218673237*Cos(5*phi) - - 0.3505558966764394*Cos(7*phi) + - 0.06575271373888199*Cos(9*phi) - - 0.0056359468919041705*Cos(11*phi) + - (5.771209617309871*Log(Cos(phi/2.)) - - 5.771209617309871*Log(Sin(phi/2.)))*Sin(phi)**12) qost(8,6) = Sin(phi)**6*(9.696916641259447*Cos(phi) + ! n=8 m=6 - 2.193094229287059*Cos(3*phi) - - 0.8592551008864449*Cos(5*phi) + - 0.19839488794352628*Cos(7*phi) - - 0.020523609097606168*Cos(9*phi) + - (-5.254043928987179*Log(Cos(phi/2.)) + - 5.254043928987179*Log(Sin(phi/2.)))*Sin(phi)**10) qost(8,8) = Sin(phi)**8*(3.0053970494755093*Cos(phi) - ! n=8 m=8 - 1.594591230881425*Cos(3*phi) + - 0.47879371162236*Cos(5*phi) - - 0.06245135368987304*Cos(7*phi) + - (3.9968866361518747*Log(Cos(phi/2.)) - - 3.9968866361518747*Log(Sin(phi/2.)))*Sin(phi)**8) qost(8,:) = qost(8,:)/st qost(9,1) = -4.107998075232154*Cos(phi)*Sqrt(Sin(phi)**2)* ! n=9 m=1 - (-0.9205149938288214 + Sin(phi)**2)* - (-0.7023117672441235 + Sin(phi)**2)* - (-0.4079422861348879 + Sin(phi)**2)* - (-0.13828286835463852 + Sin(phi)**2)* - (2.514646398761437 - 0.3827387765191469*Sin(phi)**2 + - Sin(phi)**4)*(3.2283756672875428 + 3.0517906920816182* - Sin(phi)**2 + Sin(phi)**4) qost(9,3) = 10.21134162258007*Cos(phi)*(Sin(phi)**2)**1.5* ! n=9 m=3 - (-0.9146728071791576 + Sin(phi)**2)* - (-0.6790682701752889 + Sin(phi)**2)* - (-0.3570586812035815 + Sin(phi)**2)* - (2.282752045555959 - 0.4129126943697554*Sin(phi)**2 + - Sin(phi)**4)*(2.8655463883440673 + 2.8637124529277833* - Sin(phi)**2 + Sin(phi)**4) qost(9,5) =-11.716374902027903*Cos(phi)*Sin(phi)**4 ! n=9 m=5 - *Sqrt(Sin(phi)**2)*(-0.8997067998640328 + Sin(phi)**2)* - (-0.6148470544344792 + Sin(phi)**2)* - (1.8450437598325922 - 0.4716037293271186*Sin(phi)**2 + - Sin(phi)**4)*(2.200405657268815 + 2.4861575836256304* - Sin(phi)**2 + Sin(phi)**4) qost(9,7) = 9.53629430601369*Cos(phi)*Sin(phi)**6* ! n=9 m=7 - Sqrt(Sin(phi)**2)*(-0.86169405438674 + Sin(phi)**2)* - (1.25076750372445 - 0.5516495672762026*Sin(phi)**2 + - Sin(phi)**4)*(1.353091347775423 + 1.9133436216629427* - Sin(phi)**2 + Sin(phi)**4) qost(9,9) = -0.000052377686964656096*Sqrt(Sin(phi)**2)* ! n=9 m=9 - (9799.999999999925*Cos(phi) - 21952.000000000022*Cos(3*phi) + - 20383.999999999985*Cos(5*phi) - 11872.000000000004* - Cos(7*phi) + 4760.*Cos(9*phi) - 1359.9999999999998* - Cos(11*phi) + 272.*Cos(13*phi) - 34.*Cos(15*phi) + 2.* - Cos(17*phi) + (0,0.)*Sin(phi) - (0,7.275957614183426e-12)* - Sin(3*phi) - (0,7.275957614183426e-12)*Sin(5*phi) + (0,0.)* - Sin(7*phi) +(0,0.)*Sin(9*phi) + (0,0.)*Sin(11*phi) + (0,0.)* - Sin(13*phi) + (0,0.)*Sin(15*phi) + (0,0.)*Sin(17*phi)) qost(9,:) = -qost(9,:)/st if(iodd_qnm.eq.1) then QOST(2,1) = 0.9128709291752769*Sin(phi)**4 - 0.912870929175277* | Sqrt(Sin(phi)**2) qost(2,1) = qost(2,1)/st QOST(3,2) = 1.8114220932736798*(-1 + Csc**4)*Sin(phi)**6 qost(3,2) = -qost(3,2)/st QOST(4,1) = -1.3895565288748695*Sin(phi)**8 + | Sqrt(Sin(phi)**2)*(-0.95831484749991 + | 2.3478713763747794*Sin(phi)**2) qost(4,1) = qost(4,1)/st QOST(4,3) = 2.662235902394827*Sin(phi)**8 - | 2.6622359023948277*(Sin(phi)**2)**1.5 qost(4,3) = qost(4,3)/st QOST(6,1) = 2.1052738319087974*Sin(phi)**12 - | 6.884814423269311*(-0.7645855011372953 + Sin(phi))* | (-0.5137060053527499 + Sin(phi))*Sqrt(Sin(phi)**2)* | (0.5137060053527499 + Sin(phi))*(0.7645855011372953 + | Sin(phi)) qost(6,1) = qost(6,1)/st QOST(6,3) = -4.728196355251072*Sin(phi)**12 + | 10.885847422554793*(-0.7521014330903549 + Sin(phi))* | (Sin(phi)**2)**1.5*(0.7521014330903549 + Sin(phi)) qost(6,3) = qost(6,3)/st QOST(6,5) = 4.237310362381174*Sin(phi)**12 - | 4.237310362381175*(Sin(phi)**2)**2.5 qost(6,5) = qost(6,5)/st QOST(7,6) = 4.974329856325495*(-1 + Csc**8)*Sin(phi)**14 qost(7,6) = -qost(7,6)/st QOST(8,1) = -3.223376371998216*Sin(phi)**16 + | 21.71422190007495*(-0.8051162082594958 + Sin(phi))* | (-0.7019450872852588 + Sin(phi))*(-0.41006328136624415 + | Sin(phi))*Sqrt(Sin(phi)**2)*(0.41006328136624415 +Sin(phi))* | (0.7019450872852588 + Sin(phi))*(0.8051162082594958 + | Sin(phi)) qost(8,1) = qost(8,1)/st QOST(8,3) = 7.773361906919686*Sin(phi)**16 - | 40.252589930260996*(-0.8134070954804088 + Sin(phi))* | (-0.6551930920966856 + Sin(phi))*(Sin(phi)**2)**1.5* | (0.6551930920966856 + Sin(phi))*(0.8134070954804088 + | Sin(phi)) qost(8,3) = qost(8,3)/st QOST(8,5) = Sin(phi)**4*(-8.298243772511006*Sin(phi)**12 + | Sqrt(Sin(phi)**2)*(-15.72298820054717 + | 24.02123197305818*Sin(phi)**2)) qost(8,5) = qost(8,5)/st QOST(8,7) =5.684460993638224*Sin(phi)**16 - | 5.684460993638225*(Sin(phi)**2)**3.5 qost(8,7) = qost(8,7)/st endif ! return end subroutine calqnm !----------------------------------------------------------------------- end module coef_module