/[PAMELA software]/DarthVader/OrbitalInfo/src/igrf_sub.for
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/igrf_sub.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.7 by mocchiut, Tue May 15 15:00:51 2012 UTC revision 1.11 by mocchiut, Fri Oct 31 14:02:53 2014 UTC
# Line 1  Line 1 
1          subroutine igrf_sub(xlat,xlong,year,height,  c        subroutine igrf_sub(xlat,xlong,year,height,
2       &          xl,icode,dip,dec)  c     &          xl,icode,dip,dec)
3  c----------------------------------------------------------------  c----------------------------------------------------------------
4  c   INPUT:  c   INPUT:
5  c     xlat     geodatic latitude in degrees  c     xlat     geodatic latitude in degrees
# Line 14  c          =3  an approximation is used Line 14  c          =3  an approximation is used
14  c     dip     geomagnetic inclination in degrees  c     dip     geomagnetic inclination in degrees
15  c     dec     geomagnetic declination in degress  c     dec     geomagnetic declination in degress
16  c----------------------------------------------------------------  c----------------------------------------------------------------
17    c
18        REAL              LATI,LONGI  c      REAL              LATI,LONGI
19        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  c      COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
20        SAVE /GENER/  c      SAVE /GENER/
21  C  C
22        CALL INITIZE  c      CALL INITIZE
23        ibbb=0  c      ibbb=0
24        ALOG2=ALOG(2.)  c      ALOG2=ALOG(2.)
25        ISTART=1  c      ISTART=1
26          lati=xlat  c        lati=xlat
27          longi=xlong  c        longi=xlong
28  c  c
29  C----------------CALCULATE PROFILES-----------------------------------  C----------------CALCULATE PROFILES-----------------------------------
30  c  c
31          CALL FELDCOF(YEAR,DIMO)  c        CALL FELDCOF(YEAR,DIMO)
32          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  c        CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
33          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  c        CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)
34          DIP=ASIN(BDOWN/BABS)/UMR  c        DIP=ASIN(BDOWN/BABS)/UMR
35        DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  c      DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
36        RETURN  c      RETURN
37        END  c      END
38  c  c
39  c  c
40  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
# Line 196  C--------------------------------------- Line 196  C---------------------------------------
196        SAVE /FIDB0/        SAVE /FIDB0/
197        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
198        SAVE /GENER/        SAVE /GENER/
199          REAL FLS
200  C  C
201  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
202  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
# Line 203  C-- STEQ IS STEP SIZE FOR INTEGRATION Line 204  C-- STEQ IS STEP SIZE FOR INTEGRATION
204  C  C
205        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
206        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
207          BEQU=1.E10        BEQU=1.E10
208          FLS = FL
209  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
210        RLAT=GLAT*UMR        RLAT=GLAT*UMR
211        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 306  C*****INNER LOOP (FOR QUADRATURE)      
306        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
307        ZQ=Z*Z        ZQ=Z*Z
308        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
309          IF(R.NE.R)THEN
310             FL = FLS
311             RETURN
312          ENDIF
313        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
314        RQ=R*R        RQ=R*R
315        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
# Line 677  C--   DETERMINE MAGNETIC DIPOL MOMENT AN Line 683  C--   DETERMINE MAGNETIC DIPOL MOMENT AN
683           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
684           F0 = F0 + F * F           F0 = F0 + F * F
685   1234 CONTINUE   1234 CONTINUE
686        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
687                
688        GH1(1) =  0.0        GH1(1) =  0.0
689        I=2                  I=2          
# Line 693  c      print *, "quq" Line 699  c      print *, "quq"
699           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
700           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
701           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
702           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
703           I = I+1                                                   I = I+1                                        
704           DO 9 M=1,N                                               DO 9 M=1,N                                    
705              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
706              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))
707              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
708              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
709              I=I+2              I=I+2
710   9       CONTINUE                                             9       CONTINUE                                          
711        RETURN        RETURN
# Line 947  C ERA, EREQU and ERPOL as recommended by Line 953  C ERA, EREQU and ERPOL as recommended by
953  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
954  C-----------------------------------------------------------------  C-----------------------------------------------------------------
955          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
956          CHARACTER (len=258) TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2,TP3
957          INTEGER L1,L2,L3          INTEGER L1,L2,L3
958          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
959          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23