/[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.3 by mocchiut, Tue Aug 4 14:01:27 2009 UTC revision 1.9 by mocchiut, Fri Jan 17 09:54:31 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
6  c       xlong   geodatic longitude in degrees  c     xlong    geodatic longitude in degrees
7  c       year    decimal year (year+month/12.0-0.5 or year+day-of-year/365  c     year     decimal year (year+month/12.0-0.5 or year+day-of-year/365
8  c               or 366 if leap year)  c          or 366 if leap year)
9  c       height  height in km  c     height     height in km
10  c   OUTPUT:  c   OUTPUT:
11  c       xl      L value  c     xl     L value
12  c       icode   =1  L is correct; =2  L is not correct;  c     icode     =1  L is correct; =2  L is not correct;
13  c               =3  an approximation is used  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 191  C                          APPROXIMATION Line 191  C                          APPROXIMATION
191  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
192  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
193        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
194        COMMON            X(3),H(144)        COMMON            X(3),H(196)
195        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
196        SAVE /FIDB0/        SAVE /FIDB0/
197        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 342  C-- The minimal allowable value of FI wa Line 342  C-- The minimal allowable value of FI wa
342  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
343  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
344  C  C
345  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
346  C  C
347  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
348  C  C
349  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
350  C  C
351        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
352          arg1=alog(FI)        arg1=alog(FI)
353          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
354  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
355  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
356        XX=3*arg1-arg2        XX=3*arg1-arg2
357        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
358        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
# Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN
394  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
395  C*******************************************************************  C*******************************************************************
396        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
397        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
398  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
399        ZM=P(3)                                                                  ZM=P(3)                                                          
400        FLI=P(1)*P(1)+P(2)*P(2)+1E-15        FLI=P(1)*P(1)+P(2)*P(2)+1E-15
# Line 476  C                 TO THE LOCAL GEODETIC Line 476  C                 TO THE LOCAL GEODETIC
476  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
477  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
478  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
479        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
480        CHARACTER*258      NAME        CHARACTER*258      NAME
481        INTEGER NMAX        INTEGER NMAX
482        REAL TIME        REAL TIME
483        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
484        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
485        SAVE/MODEL/        SAVE/MODEL/
486        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 581  C--------------------------------------- Line 581  C---------------------------------------
581        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
582        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
583  C     ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
584  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
585        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
586        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
587        INTEGER L1,L2,L3        INTEGER L1,L2,L3
588        INTEGER NMAX        INTEGER NMAX
# Line 615  c      print *, "qua" Line 615  c      print *, "qua"
615  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
616  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
617  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
618  c     WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
619  c     WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
620  c     WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
621  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
622        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
623  c      c    
624  c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
625  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
# Line 652  c     L = (IYEA - 1945)/5 + 1 Line 652  c     L = (IYEA - 1945)/5 + 1
652        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
653        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
654        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
655    c      WRITE(*,*) FIL1
656    c      WRITE(*,*) FIL2
657  c      print *, "que"  c      print *, "que"
658  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
659        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
660        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
661    c      print *, "quessss"
662        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
663        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
664  c      print *, "quj"  c      print *, "quj"
# Line 674  C--   DETERMINE MAGNETIC DIPOL MOMENT AN Line 677  C--   DETERMINE MAGNETIC DIPOL MOMENT AN
677           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
678           F0 = F0 + F * F           F0 = F0 + F * F
679   1234 CONTINUE   1234 CONTINUE
680        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
681                
682        GH1(1) =  0.0        GH1(1) =  0.0
683        I=2                  I=2          
# Line 690  c      print *, "quq" Line 693  c      print *, "quq"
693           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
694           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
695           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
696           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
697           I = I+1                                                   I = I+1                                        
698           DO 9 M=1,N                                               DO 9 M=1,N                                    
699              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
700              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))
701              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
702              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
703              I=I+2              I=I+2
704   9       CONTINUE                                             9       CONTINUE                                          
705        RETURN        RETURN
# Line 944  C ERA, EREQU and ERPOL as recommended by Line 947  C ERA, EREQU and ERPOL as recommended by
947  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
948  C-----------------------------------------------------------------  C-----------------------------------------------------------------
949          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
950          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=258) TP1,TP2,TP3
951          INTEGER L1,L2,L3          INTEGER L1,L2,L3
952          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
953          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23