/[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.8 by mocchiut, Thu Jan 16 15:29:26 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 55  C*************************************** Line 55  C***************************************
55  C  C
56  C  C
57        SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)        SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)
58          IMPLICIT REAL(8)(A-H)
59          IMPLICIT REAL(8)(O-Z)
60  C--------------------------------------------------------------------  C--------------------------------------------------------------------
61  C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE  C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE
62  C  C
# Line 70  C          BEQU   MAGNETIC FIELD STRENGT Line 72  C          BEQU   MAGNETIC FIELD STRENGT
72  C          RR0    EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS  C          RR0    EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS
73  C          BDEL   FINAL ACHIEVED ACCURACY  C          BDEL   FINAL ACHIEVED ACCURACY
74  C--------------------------------------------------------------------  C--------------------------------------------------------------------
75          REAL(8) P
76        DIMENSION         P(8,4),SP(3)        DIMENSION         P(8,4),SP(3)
77        LOGICAL           VALUE        LOGICAL           VALUE
78        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
# Line 154  C        DO 1111 I=1,8 Line 157  C        DO 1111 I=1,8
157  C  C
158  C  C
159        SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)        SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)
160          IMPLICIT REAL(8)(A-H)
161          IMPLICIT REAL(8)(O-Z)
162  C--------------------------------------------------------------------  C--------------------------------------------------------------------
163  C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE  C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE
164  C AND GEMAGNETIC FIELD MODEL.  C AND GEMAGNETIC FIELD MODEL.
# Line 190  C                          WHICH ACCURAT Line 195  C                          WHICH ACCURAT
195  C                          APPROXIMATION IS USED.  C                          APPROXIMATION IS USED.
196  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
197  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
198          REAL(8) AQUAD,BQUAD,ERA
199          REAL(8) CT,ST,D,ALT,X,RQ,P,U,STEP
200        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
201        COMMON            X(3),H(144)        COMMON            X(3),H(196)
202        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
203        SAVE /FIDB0/        SAVE /FIDB0/
204        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     ERA,AQUAD,BQUAD,UMR
205        SAVE /GENER/        SAVE /GENER/
206  C  C
207  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
# Line 349  C Line 356  C
356  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
357  C  C
358        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
359        arg1=alog(FI)        arg1=dlog(FI)
360        arg2=alog(DIMOB0)        arg2=dlog(DIMOB0)
361  c    arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
362  c    if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
363        XX=3*arg1-arg2        XX=3*arg1-arg2
# Line 378  c  771 GG=3.33338E-1*XX+3.0062102E-1     Line 385  c  771 GG=3.33338E-1*XX+3.0062102E-1    
385       1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0                   1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0            
386        GOTO777                                                                  GOTO777                                                          
387    776 GG=XX-3.0460681E0                                                    776 GG=XX-3.0460681E0                                                
388    777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0)    777 FL=EXP(dLOG((1.+EXP(GG))*DIMOB0)/3.0)
389        RETURN                                                                    RETURN                                                            
390  C*****APPROXIMATION FOR HIGH VALUES OF L.                                C*****APPROXIMATION FOR HIGH VALUES OF L.                              
391  30    ICODE=3                                                            30    ICODE=3                                                          
# Line 389  C*****APPROXIMATION FOR HIGH VALUES OF L Line 396  C*****APPROXIMATION FOR HIGH VALUES OF L
396  C  C
397  C  C
398        SUBROUTINE STOER(P,BQ,R)                                                  SUBROUTINE STOER(P,BQ,R)                                          
399          IMPLICIT REAL(8)(A-H)
400          IMPLICIT REAL(8)(O-Z)
401    
402  C*******************************************************************  C*******************************************************************
403  C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG                *  C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG                *
404  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
405  C*******************************************************************  C*******************************************************************
406          REAL(8) P,ZM,FLI,WR,XM,R,YM,XI,DR
407        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
408        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
409  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
410        ZM=P(3)                                                                  ZM=P(3)                                                          
411        FLI=P(1)*P(1)+P(2)*P(2)+1E-15        FLI=P(1)*P(1)+P(2)*P(2)+1E-15
# Line 434  C*****FORM SLOWLY VARYING EXPRESSIONS   Line 445  C*****FORM SLOWLY VARYING EXPRESSIONS  
445  C  C
446  C  C
447        SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)                  SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)          
448          IMPLICIT REAL(8)(A-H)
449          IMPLICIT REAL(8)(O-Z)
450  C-------------------------------------------------------------------  C-------------------------------------------------------------------
451  C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL  C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
452  C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,  C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
# Line 476  C                 TO THE LOCAL GEODETIC Line 489  C                 TO THE LOCAL GEODETIC
489  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
490  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
491  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
492        DIMENSION         V(3),B(3),G(144)        REAL(8) G,H,X,Z,F,S,T,XI,Y,XXX,BXXX,BYYY,BZZZ,B
493          REAL(8) YYY,ZZZ,BABS,BEAST,BRHO,BNORTH,BDOWN,ST,CT,CP,SP
494          REAL(8) ERA, AQUAD, BQUAD,ALT,RHO,RQ,D
495          DIMENSION         V(3),B(3),G(196)
496        CHARACTER*258      NAME        CHARACTER*258      NAME
497        INTEGER NMAX        INTEGER NMAX
498        REAL TIME        REAL TIME
499        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
500          
501        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
502        SAVE/MODEL/        SAVE/MODEL/
503        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     ERA,AQUAD,BQUAD,UMR
504        SAVE/GENER/        SAVE/GENER/
505  C  C
506  C-- IS RECORDS ENTRY POINT  C-- IS RECORDS ENTRY POINT
# Line 566  c      print *,' I ',I Line 583  c      print *,' I ',I
583  C  C
584  C  C
585        SUBROUTINE FELDCOF(YEAR,DIMO)        SUBROUTINE FELDCOF(YEAR,DIMO)
586          IMPLICIT REAL(8)(A-H)
587          IMPLICIT REAL(8)(O-Z)
588  C------------------------------------------------------------------------  C------------------------------------------------------------------------
589  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
590  C      C    
# Line 581  C--------------------------------------- Line 600  C---------------------------------------
600        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
601        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
602  C     ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
603  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
604        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        REAL(8) GH1, GH2, GHA
605          DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
606        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
607          DOUBLE PRECISION DIMO
608        INTEGER L1,L2,L3        INTEGER L1,L2,L3
609        INTEGER NMAX        INTEGER NMAX
610          REAL YEAR
611        REAL TIME        REAL TIME
612        CHARACTER *258 P1,P2,P3        CHARACTER *258 P1,P2,P3
613          REAL(8) AQUAD,BQUAD,ERAD
614        COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/PPATH/ L1,L2,L3,P1, P2, P3
615        SAVE/PPATH/        SAVE/PPATH/
616        COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
617        SAVE/MODEL/        SAVE/MODEL/
618        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/GENER/   ERAD,AQUAD,BQUAD,UMR
619        SAVE/GENER/        SAVE/GENER/
620  C     ### updated to 2005  C     ### updated to 2005
621  C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80
# Line 615  c      print *, "qua" Line 638  c      print *, "qua"
638  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
639  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
640  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
641  c     WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
642  c     WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
643  c     WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
644  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
645        DATA   DTEMOD / 2005., 2010., 2015./        DATA   DTEMOD / 2005., 2010., 2015./
646  c      c    
# Line 652  c     L = (IYEA - 1945)/5 + 1 Line 675  c     L = (IYEA - 1945)/5 + 1
675        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
676        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
677        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
678    c      WRITE(*,*) FIL1
679    c      WRITE(*,*) FIL2
680  c      print *, "que"  c      print *, "que"
681  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
682        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
683        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
684    c      print *, "quessss"
685        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
686        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
687  c      print *, "quj"  c      print *, "quj"
# Line 704  c      print *, "quq" Line 730  c      print *, "quq"
730  C      C    
731  C      C    
732        SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
733                                                                                          IMPLICIT REAL(8)(A-H)
734          IMPLICIT REAL(8)(O-Z)
735  C ===============================================================                C ===============================================================              
736  C                                                                                C                                                                              
737  C     Version 1.01                                                  C     Version 1.01                                                
# Line 733  C     Line 760  C    
760  C     ===============================================================                C     ===============================================================              
761                
762        CHARACTER  FSPEC*(*), FOUT*258        CHARACTER  FSPEC*(*), FOUT*258
763          REAL(8) GH,ERAD
764        DIMENSION       GH(*)                                                DIMENSION       GH(*)                                        
765  C     ---------------------------------------------------------------                C     ---------------------------------------------------------------              
766  C     Open coefficient file. Read past first header record.          C     Open coefficient file. Read past first header record.        
# Line 796  C Line 824  C
824  C  C
825          SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,                    SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,          
826       1                        NMAX2, GH2, NMAX, GH)                         1                        NMAX2, GH2, NMAX, GH)                  
827                                                                                            IMPLICIT REAL(8)(A-H)
828            IMPLICIT REAL(8)(O-Z)
829            REAL DATE
830  C ===============================================================                C ===============================================================              
831  C                                                                                C                                                                              
832  C       Version 1.01                                                  C       Version 1.01                                                
# Line 824  C       USGS, MS 964, Box 25046 Federal Line 854  C       USGS, MS 964, Box 25046 Federal
854  C                                                                                C                                                                              
855  C ===============================================================                C ===============================================================              
856                                                                                                                                                                    
857            REAL(8) GH1, GH2, GH
858          DIMENSION       GH1(*), GH2(*), GH(*)                                  DIMENSION       GH1(*), GH2(*), GH(*)                        
859                                                                                                                                                                    
860  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
# Line 863  C Line 894  C
894  C  C
895          SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,                    SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,          
896       1                        GH2, NMAX, GH)                                 1                        GH2, NMAX, GH)                          
897            IMPLICIT REAL(8)(A-H)
898            IMPLICIT REAL(8)(O-Z)
899            REAL DATE
900                                                                                                                                                                    
901  C ===============================================================                C ===============================================================              
902  C                                                                                C                                                                              
# Line 891  C       USGS, MS 964, Box 25046 Federal Line 925  C       USGS, MS 964, Box 25046 Federal
925  C                                                                                C                                                                              
926  C ===============================================================                C ===============================================================              
927                                                                                                                                                                    
928            REAL(8) GH1, GH2, GH
929          DIMENSION       GH1(*), GH2(*), GH(*)                                    DIMENSION       GH1(*), GH2(*), GH(*)                          
930                                                                                                                                                                    
931  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
# Line 929  C -------------------------------------- Line 964  C --------------------------------------
964  C  C
965  C  C
966          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)
967            IMPLICIT REAL(8)(A-H)
968            IMPLICIT REAL(8)(O-Z)
969  C----------------------------------------------------------------  C----------------------------------------------------------------
970  C Initializes the parameters in COMMON/GENER/  C Initializes the parameters in COMMON/GENER/
971  C  C
# Line 944  C ERA, EREQU and ERPOL as recommended by Line 981  C ERA, EREQU and ERPOL as recommended by
981  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
982  C-----------------------------------------------------------------  C-----------------------------------------------------------------
983          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
984          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=258) TP1,TP2,TP3
985          INTEGER L1,L2,L3          INTEGER L1,L2,L3
986          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
987            REAL(8) AQUAD,BQUAD,ERA
988    
989          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
990          SAVE/PPATH/          SAVE/PPATH/
991    
992          COMMON/GENER/UMR,ERA,AQUAD,BQUAD          COMMON/GENER/ERA,AQUAD,BQUAD,UMR
993          SAVE/GENER/          SAVE/GENER/
994    
995          L1 = TL1          L1 = TL1
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23