/[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.5 by mocchiut, Wed Nov 23 21:19:34 2011 UTC revision 1.12 by mocchiut, Mon Jan 19 12:32:13 2015 UTC
# Line 1  Line 1 
         subroutine igrf_sub(xlat,xlong,year,height,  
      &          xl,icode,dip,dec)  
 c----------------------------------------------------------------  
 c   INPUT:  
 c     xlat     geodatic latitude in degrees  
 c     xlong    geodatic longitude in degrees  
 c     year     decimal year (year+month/12.0-0.5 or year+day-of-year/365  
 c          or 366 if leap year)  
 c     height     height in km  
 c   OUTPUT:  
 c     xl     L value  
 c     icode     =1  L is correct; =2  L is not correct;  
 c          =3  an approximation is used  
 c     dip     geomagnetic inclination in degrees  
 c     dec     geomagnetic declination in degress  
 c----------------------------------------------------------------  
   
       REAL              LATI,LONGI  
       COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  
       SAVE /GENER/  
1  C  C
       CALL INITIZE  
       ibbb=0  
       ALOG2=ALOG(2.)  
       ISTART=1  
         lati=xlat  
         longi=xlong  
 c  
 C----------------CALCULATE PROFILES-----------------------------------  
 c  
         CALL FELDCOF(YEAR,DIMO)  
         CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  
         CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  
         DIP=ASIN(BDOWN/BABS)/UMR  
       DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  
       RETURN  
       END  
 c  
 c  
2  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
3  C  C
4  C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2    C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2  
# Line 191  C                          APPROXIMATION Line 153  C                          APPROXIMATION
153  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
154  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
155        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
156        COMMON            X(3),H(144)        COMMON            X(3),H(196)
157        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
158        SAVE /FIDB0/        SAVE /FIDB0/
159        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
160        SAVE /GENER/        SAVE /GENER/
161          REAL FLS
162  C  C
163  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
164  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 166  C-- STEQ IS STEP SIZE FOR INTEGRATION
166  C  C
167        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
168        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
169          BEQU=1.E10        BEQU=1.E10
170          FLS = FL
171  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
172        RLAT=GLAT*UMR        RLAT=GLAT*UMR
173        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 268  C*****INNER LOOP (FOR QUADRATURE)      
268        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
269        ZQ=Z*Z        ZQ=Z*Z
270        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
271          IF(R.NE.R)THEN
272             FL = FLS
273             RETURN
274          ENDIF
275        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
276        RQ=R*R        RQ=R*R
277        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
# Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN Line 362  C* SUBROUTINE USED FOR FIELD LINE TRACIN
362  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
363  C*******************************************************************  C*******************************************************************
364        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
365        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
366  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
367        ZM=P(3)                                                                  ZM=P(3)                                                          
368        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 444  C                 TO THE LOCAL GEODETIC
444  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
445  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
446  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
447        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
448        CHARACTER*258      NAME        CHARACTER*258      NAME
449        INTEGER NMAX        INTEGER NMAX
450        REAL TIME        REAL TIME
451        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
452        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
453        SAVE/MODEL/        SAVE/MODEL/
454        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 580  C     ### updated to IGRF-2005 version - Line 548  C     ### updated to IGRF-2005 version -
548  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
549        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
550        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
551  C     ### FILMOD, DTEMOD arrays +1        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(2)
 c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  
       DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)  
552        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
553        INTEGER L1,L2,L3        INTEGER L1,L2,I1
554        INTEGER NMAX        INTEGER NMAX
555        REAL TIME        REAL TIME
556        CHARACTER *258 P1,P2,P3        CHARACTER *258 P1,P2
557        COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/PPATH/ I1,L1,L2,P1,P2
558        SAVE/PPATH/        SAVE/PPATH/
559        COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
560        SAVE/MODEL/        SAVE/MODEL/
561        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
562        SAVE/GENER/        SAVE/GENER/
 C     ### updated to 2005  
 C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  
         
 c     COEFPATH  = 'OrbitalInfo/src/'  
 c     COEF1 = 'dgrf00.dat'  
 c     COEF2 = 'igrf05.dat'  
 c     COEF3 = 'igrf05s.dat'  
 c     COEF1 = COEFPATH(1:16)//COEF1  
 c     COEF2 = COEFPATH(1:16)//COEF2  
 c     COEF3 = COEFPATH(1:16)//COEF3  
 c     FILMOD(1) = COEF1  
 c     FILMOD(2) = COEF2  
 c     FILMOD(3) = COEF3  
563  c      print *, "qui"  c      print *, "qui"
564        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
565        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
       FILMOD(3) = P3(1:L3)  
566  c      print *, "qua"  c      print *, "qua"
 c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  
 c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  
 c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  
 c     WRITE(*,*) FILMOD(1)  
 c     WRITE(*,*) FILMOD(2)  
 c     WRITE(*,*) FILMOD(3)  
 c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  
       DATA   DTEMOD / 2005., 2010., 2015./  
 c      
 c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              
 c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        
 c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',        
 c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',  
 c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/  
 c     DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,  
 c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        
 C      
 C     ### numye = numye + 1 ; is number of years represented by IGRF  
 C      
 c     NUMYE=13  
       NUMYE=2  
 c      print *, "quo"  
         
567  C      C    
568  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
569  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
# Line 644  C     Line 573  C    
573  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
574        TIME = YEAR        TIME = YEAR
575        IYEA = INT(YEAR/5.)*5        IYEA = INT(YEAR/5.)*5
576  c     L = (IYEA - 1945)/5 + 1        L = IYEA + 5
577        L = (IYEA - 2000)/5 + 1  C
578        IF(L.LT.1) L=1        DTE1 = REAL(IYEA)  
579        IF(L.GT.NUMYE) L=NUMYE                FIL1 = FILMOD(1)  
580        DTE1 = DTEMOD(L)          DTE2 = REAL(L)
581        FIL1 = FILMOD(L)          FIL2 = FILMOD(2)
582        DTE2 = DTEMOD(L+1)  c      print *,'IYEA ',IYEA,' L ',L
583        FIL2 = FILMOD(L+1)  c      WRITE(*,*) FIL1
584    c      WRITE(*,*) FIL2
585  c      print *, "que"  c      print *, "que"
586  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
587        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
588        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
589    c      print *, "quessss"
590        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
591        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
592  c      print *, "quj"  c      print *, "quj"
593  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
594        IF (L .LE. NUMYE-1) THEN                                IF (I1.EQ.0) THEN      
595           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
596       1        NMAX2, GH2, NMAX, GHA)                               1        NMAX2, GH2, NMAX, GHA)                        
597        ELSE                      ELSE              
# Line 674  C--   DETERMINE MAGNETIC DIPOL MOMENT AN Line 605  C--   DETERMINE MAGNETIC DIPOL MOMENT AN
605           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
606           F0 = F0 + F * F           F0 = F0 + F * F
607   1234 CONTINUE   1234 CONTINUE
608        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
609                
610        GH1(1) =  0.0        GH1(1) =  0.0
611        I=2                  I=2          
# Line 690  c      print *, "quq" Line 621  c      print *, "quq"
621           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
622           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
623           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
624           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
625           I = I+1                                                   I = I+1                                        
626           DO 9 M=1,N                                               DO 9 M=1,N                                    
627              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
628              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))
629              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
630              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
631              I=I+2              I=I+2
632   9       CONTINUE                                             9       CONTINUE                                          
633        RETURN        RETURN
# Line 928  C -------------------------------------- Line 859  C --------------------------------------
859          END                                                                      END                                                            
860  C  C
861  C  C
862          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)          SUBROUTINE INITIZE(ISSEC,TP1,TL1,TP2,TL2)
863  C----------------------------------------------------------------  C----------------------------------------------------------------
864  C Initializes the parameters in COMMON/GENER/  C Initializes the parameters in COMMON/GENER/
865  C  C
# Line 943  C Line 874  C
874  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
875  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
876  C-----------------------------------------------------------------  C-----------------------------------------------------------------
877          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,ISSEC
878          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2
879          INTEGER L1,L2,L3          INTEGER L1,L2
880          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2
881          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ I1,L1,L2,P1,P2
882          SAVE/PPATH/          SAVE/PPATH/
883    
884          COMMON/GENER/UMR,ERA,AQUAD,BQUAD          COMMON/GENER/UMR,ERA,AQUAD,BQUAD
885          SAVE/GENER/          SAVE/GENER/
886    
887            I1 = ISSEC
888          L1 = TL1          L1 = TL1
889          L2 = TL2          L2 = TL2
890          L3 = TL3  
891          P1 = TP1(1:L1)          P1 = TP1(1:L1)
892          P2 = TP2(1:L2)          P2 = TP2(1:L2)
         P3 = TP3(1:L3)  
893    
894          ERA=6371.2          ERA=6371.2
895          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23