/[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.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 196  C--------------------------------------- Line 158  C---------------------------------------
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 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(196),GH2(196),GHA(196),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
       FIL2 = FILMOD(L+1)  
583  c      WRITE(*,*) FIL1  c      WRITE(*,*) FIL1
584  c      WRITE(*,*) FIL2  c      WRITE(*,*) FIL2
585  c      print *, "que"  c      print *, "que"
# Line 663  c      print *, "quessss" Line 591  c      print *, "quessss"
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 677  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 693  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 931  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 946  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 (len=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.7  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23