/[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.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
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
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 342  C-- The minimal allowable value of FI wa Line 348  C-- The minimal allowable value of FI wa
348  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.
349  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
350  C  C
351  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
352  C  C
353  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
354  C  C
355  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
356  C  C
357        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
358          arg1=alog(FI)        arg1=alog(FI)
359          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
360  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
361  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
362        XX=3*arg1-arg2        XX=3*arg1-arg2
363        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
364        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
# Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN Line 400  C* SUBROUTINE USED FOR FIELD LINE TRACIN
400  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
401  C*******************************************************************  C*******************************************************************
402        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
403        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
404  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
405        ZM=P(3)                                                                  ZM=P(3)                                                          
406        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 482  C                 TO THE LOCAL GEODETIC
482  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
483  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
484  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
485        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
486        CHARACTER*258      NAME        CHARACTER*258      NAME
487        INTEGER NMAX        INTEGER NMAX
488        REAL TIME        REAL TIME
489        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
490        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
491        SAVE/MODEL/        SAVE/MODEL/
492        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 581  C--------------------------------------- Line 587  C---------------------------------------
587        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
588        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
589  C     ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
590  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
591        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
592        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
593        INTEGER L1,L2,L3        INTEGER L1,L2,L3
594        INTEGER NMAX        INTEGER NMAX
# Line 615  c      print *, "qua" Line 621  c      print *, "qua"
621  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
622  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
623  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
624  c     WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
625  c     WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
626  c     WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
627  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
628        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
629  c      c    
630  c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
631  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 658  c     L = (IYEA - 1945)/5 + 1
658        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
659        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
660        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
661    c      WRITE(*,*) FIL1
662    c      WRITE(*,*) FIL2
663  c      print *, "que"  c      print *, "que"
664  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
665        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
666        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
667    c      print *, "quessss"
668        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
669        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
670  c      print *, "quj"  c      print *, "quj"
# Line 674  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 690  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 944  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 *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.3  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23