program Sun_fisheye implicit none c c This program combines the anti-sun centered ecliptic sky maps made by c program gegenscheinmovies.for, into an anti-sun centered ecliptic fisheye c This is a square 180 x 180 degrees c character*58 filename1,filename2,blankk character*62 filename3,blank3 character*46 outfile,mapfile character*10 head(5) character*12 XY character*3 DOY character*2 cs character*1 i100,i10,i1 integer*4 i,imax,j,k,jra,kdec,icount,imap,idoy real*4 frame1(360,360),frame2(720,360),x,y,r,phi,rLng,rLat,xsq,cmap(360,360) real*4 amap(360,360),arg1,arg2,argsum1,argsum2,corr32,corr3,eps!,arg10,arg20 data imax /1947/ data blankk /' '/ data blank3 /' '/ data cs /'cm'/ data mapfile /'g:\GegenscheinXStars\c3avg_2003through2008.grd'/ data outfile /'g:\GegenscheinXStars\cmecl_200X_DOY_HHMMSS.grd'/ data imap /1/ !zero means make no maps c open (10,file='G:\GegenscheinXStars\filenames123.txt',readonly) open (14,file='G:\GegenscheinXStars\Cam2-3MatchFlatMap.txt') c c Bring in camera 3 correction map and roll it off c Here the main roll off is a gaussian for eps > 50 deg, plus an c aggressive longitude rolloff starting at 60 degrees to kill some c unpleasant artifacts in the map c open (12,file=mapfile,readonly) read(12,1)head do j=1,360 read(12,11)(cmap(i,j),i=1,360) enddo close(12) do j=1,360 arg2=cosd(float(j-180)/2.) do i=1,360 arg1=cosd(float(i-180)/2.) eps=acosd(arg1*arg2) if(eps.gt.50.)then eps=eps-50. cmap(i,j)=cmap(i,j)*exp(-(eps*eps)/600.) endif if(arg1.lt.0.5)then !here reduce some longitude artifacts... cmap(i,j)=cmap(i,j)*(2.*arg1)**2 endif enddo enddo c open(12,file='g:\gegenscheinxstars\c3avg_2003-8rolledoff.grd') write(12,1)head do j=1,360 write(12,11)(cmap(i,j),i=1,360) enddo close(12) c if(imax.ne.0)go to 98 !******use this when only rolled-off map is desired to be made*** c do i=1,imax if(mod(i,100).eq.0)print *,i read(10,10)filename1,filename2,filename3 10 format(6x,a58,9x,a58,9x,a62) c print *,filename1,filename2,filename3 if(filename3.eq.blank3.and.filename2.eq.blankk)go to 99 if(i.eq.1034)go to 99 !Kill DOY 155, 2006 very little sky coverage if(i.eq.1641.or.i.eq.1642)go to 99 !Kill DOYs 121,122, 2008 "hot" c if(i.ge.1710.and.i.le.1722)go to 99 !Kill DOYs 206-218, 2008 "bad patterns", recovered c c Zero Fisheye map c do j=1,360 do k=1,360 amap(j,k)=0. enddo enddo c c Read files c if(filename3.eq.blank3)then do k=1,360 do j=1,360 frame1(j,k)=0. enddo enddo write(XY(1:12),'(a12)')filename2(43:54) else open(11,file=filename3,readonly) read(11,1)head 1 format(a10) do k=1,360 read(11,11)(frame1(j,k),j=1,360) 11 format(20f8.1) do j=1,360 if(frame1(j,k).ne.0..and.frame1(j,k).lt.10000.)frame1(j,k)=frame1(j,k)-cmap(j,k) enddo enddo close(11) write(XY(1:12),'(a12)')filename3(47:58) endif write(outfile(31:42),'(a12)')XY write(DOY(1:3),'(a3)')XY(3:5) c Decoding DOY and evaluate empirical DOY-dependent offset for camera 3 write(i100,'(a1)')XY(3:3) write(i10,'(a1)')XY(4:4) write(i1,'(a1)')XY(5:5) idoy=100*(ichar(i100)-48)+10*(ichar(i10)-48)+ichar(i1)-48 xsq=float(idoy)*float(idoy) corr3 =6.7*cos(float(idoy)/55.+2.4+0.00002*xsq)+1.73-0.03*float(idoy) !flattens cam 3 alone c corr32=7.8*cos(float(idoy)/55.+2.4+0.00002*xsq)-2.73-0.03*float(idoy) !matches cams 2 & 3, before map subtraction corr32=7.8*cos(float(idoy)/55.+2.4+0.00002*xsq)-6.06-0.03*float(idoy) !matches cams 2 & 3, after map subtraction c c Include ad hoc mask-in corrections: ***************depends upon sequence number in list******************** c if(i.ge.25.and.i.le.42)corr32=corr32-18. !ad hoc correction (no mask) early in SMEI, wasn't needed before update... if(i.gt. 972.and.i.le. 996)corr32=corr32 -0.5*float(i-966) if(i.gt.1136.and.i.le.1250)corr32=corr32 +0.22*float(i-1136) if(i.gt.1250.and.i.le.1325)corr32=corr32+28.-0.04*float(i-1250) if(i.gt.1434.and.i.le.1655)then if(i.lt.1581.or.i.gt.1591)corr32=corr32+52.-0.15*float(i-1435) endif if(i.ge.1750)then if(i.le.1901.or.i.ge.1915)corr32=corr32+70. endif c if(filename2.eq.blankk)then do k=1,360 do j=1,720 frame2(j,k)=0. enddo enddo else open(11,file=filename2,readonly) read(11,1)head do k=1,360 read(11,11)(frame2(j,k),j=1,720) enddo close(11) endif c c Load fisheye map: if both cameras are present in a bin, use average (or the closer-to-zero one, commented out below) c icount=0 argsum1=0. argsum2=0. do j=1,360 y=89.75-0.5*float(j-1) do k=1,360 arg1=0. arg2=0. x=89.75-0.5*float(k-1) r=sqrt(x*x+y*y) c if(r.gt.2.)go to 100 phi=atan2d(y,x)-90. if(phi.lt.0.)phi=phi+360. rLng=phi rLat=90.-r call rotate(0.,-90.,0.,rLng,rLat) rLng=rLng+180. !antisolar: note this should be commented out for Sun-centered maps if(rLng.ge.360.)rLng=rLng-360. c jra=nint(0.25+2.*rLng) jra=jra-180 !camera 3 only kdec=nint(0.25+2.*(rLat+90.)) c print *,'jra = ',jra,kdec,x,y,phi,rLng,rLat if(jra.gt.0.and.jra.le.360)then if(frame1(jra,kdec).ne.0..and.frame1(jra,kdec).lt.10000.)then arg1=frame1(jra,kdec)-corr32 endif endif jra=jra+540 !camera 2 only if(jra.gt.720)jra=jra-720 if(jra.gt.0.and.jra.le.720)then if(frame2(jra,kdec).ne.0..and.frame2(jra,kdec).lt.10000.)arg2=frame2(jra,kdec) endif if(arg1.eq.0..and.arg2.eq.0.)go to 100 if(arg1.ne.0..and.arg2.eq.0.)then amap(k,j)=arg1 go to 100 endif if(arg1.eq.0..and.arg2.ne.0.)then amap(k,j)=arg2 go to 100 endif if(arg1.ne.0..and.arg2.ne.0.)then amap(k,j)=(arg1+arg2)/2.0 !average icount=icount+1 c if(arg1.gt.100.)print *,filename3,arg1,jra,kdec argsum1=argsum1+arg1 argsum2=argsum2+arg2 c c Chooses closest-to-zero one c c arg10=abs(arg1) c arg20=abs(arg2) c if(arg10.gt.arg20)then c amap(k,j)=arg2 c else c amap(k,j)=arg1 c endif endif 100 continue enddo enddo c c Output fisheye map c if(imap.ne.0)then print *,i,' writing ',outfile open(12,file=outfile) write(12,12) 12 format('DSAA'/'360 360'/'-90 90'/'-90 90'/'-50 50') do j=1,360 write(12,11)(amap(k,j),k=1,360) enddo close(12) endif if(icount.gt.0)then argsum1=argsum1/float(icount) argsum2=argsum2/float(icount) c write(14,14)i,DOY,argsum1,argsum2,argsum1-argsum2-corr32,argsum1-corr3,icount write(14,14)i,DOY,argsum1,argsum2,argsum1-argsum2,icount 14 format(i4,2x,a3,2x,3f10.3,i7) endif 99 enddo close(10) close(14) 98 stop 'done' end c********************************************************************************* subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT) implicit none c INPUTS: c ALFA real phase angle of the pole Z(new) in old coordinate system c BETA real polar angle of the pole Z(new) in old coordinate system c GAMMA real phase angle of X(new) in coordinate system (2) c PHI real phase angle in old coordinate system c RLAT real latitude in old coordinate system (-90<=RLAT<=90) c OUTPUTS: c PHI real phase angle in new coordinate system (0<=PHI<360) c RLAT real latitude in new coordinate system real ALFA real BETA real GAMMA real PHI real RLAT double precision SINB double precision COSB double precision SINL double precision COSL double precision SINP double precision COSP double precision COST ! Needed for g77 when using goniod.f ! What happens if these are available as intrinsic functions?? double precision dsind double precision dcosd double precision dacosd double precision datan2d if (abs(RLAT) .gt. 90.0)print *,'Bummer: RLAT = ',RLAT !call Say('rotate','E','invalid','latitude') PHI = PHI-ALFA ! Rotation (1) SINB = dsind(dble(BETA)) ! Rotation (2) COSB = dcosd(dble(BETA)) SINL = dsind(dble(RLAT)) COSL = dcosd(dble(RLAT)) SINP = dsind(dble(PHI )) COSP = dcosd(dble(PHI )) COST = SINL*COSB+COSL*SINB*COSP ! cos(polar angle) if (dabs(COST) .ne. 1.0d0) then ! changed abs --> dabs from Paul's, probably doesn't matter... SINP = COSL*SINP COSP = COSL*COSB*COSP-SINL*SINB else !Point on new Z-axis; PHI not defined !------- ! The following statement is inserted for the benefit of the ! synoptic map drawing program. ! If BETA=0 then SINB=0, COSB=+/-1, i.e. the z-axis stays the same ! or is flipped by 180 degrees. The input phase angle for points on ! the z-axis will remain the same or is changed to 180-(phase angle) ! (just as for points not on the z-axis). if (SINB .eq. 0.0d0) COSP = COSB*COSP end if PHI = datan2d(SINP,COSP) PHI = PHI-GAMMA ! Rotation (3) PHI = amod(amod(PHI,360.0)+360.0,360.0)!0<=PHI<360 RLAT = 90.0d0-dacosd(COST) ! acosd returns between 0 and 180 deg return end c*********************************************************************************