c CODE_TO_DONES/mcnewd_asyobs.for c calculation of azimuthal asymmetry 3/12/99 c 25.09.2000: use same subroutines as with mcnew_test c 9.01.2001: Henyey-Greenstein c 17.01.2001: EULER c 17.01.2001: Illumination by Saturn c 21.01.2001: Callisto, power-law c 26.01.2001: correct observer direction -> mcnew_asyobs.for c 13.06.2001: double precision version c 04.01.2002: added comment for Luke c mcnew_sub.inc = commons for transfer between subroutines c mcnew_main.inc = common for transfer between main/subroutines c includes variables related to Saturn's illumination c**************************************************************** program mcnewd_asyobs c**************************************************************** c driver for calculating azimuthal asymmetry: c loops over position files and azimuths c uses mcnewd_sub.for for the calculation of scattering c from a single position-file for each sun/observer direction c corresponding to the given ring longitude c output of results via allfile created in this program c use mc_asyobs_read.pro for reading the output c***************************************************************** c HS 01/12/99 - 19/12/01 implicit double precision (a-h,o-z) include 'mcnew_main.inc' parameter(nnt=1000) character*60 posfile,parfile,outfile,allfile,runide double precision testfii(nnt),testlat(nnt) double precision testw(nnt),testd(nnt),testn(nnt) double precision testm(nnt,10) character*60 epsfile,wwwfile double precision iscat_mean,npot_mean,nbox_mean double precision dlobs double precision azitable(100),suntable(100),obstable(100) double precision tablew(100),tabled(100),tablen(100) double precision tablem(100,10) double precision table_fdirect(100) character*8 date character*10 time character*5 zone integer*4 values(8) character*30 smet,smodel c posfile = posdat-file containing particle positions c parfile = parameter-file (junk -> defaults) c outfile = output of indirect+direct for a single pos/view (junk) c allfile = output for all pos/view c epsfile = output of weights + fiilat (junk) c wwwfile = output of weights (junk) c sredu = reduction factor for particle radius (sredu < 1) c zredu = reduction factor for vertical coordinates (zredu > 1) c idum = seed ar=0.017453293d0 ra=57.29577951d0 c******************************************************** c READ INPUT VALUES c******************************************************** c runide = run-identification c len = length of runide c positions to be used runide.posdatNN - runide.posdatMM c NN=pos1, MM=pos2 c BSUN,BOBS = elevation of sun, observer c DLOBS=LSUN-LOBS = longitude of sun with respect to observer c AZI1 - AZI2 azimuth values c NPHOT number of photons/view c OUTPUT: c testw,testn,testd,testm c to convert testw to I/F: c multiply by sin(bsun)/nphot/sin(bobs) read(*,*) runide,len read(*,*) ipos1,ipos2 read(*,*) bobs,bsun,dlobs npos=pos2-pos1+1 read(*,*) azi1,azi2,nazi azimin=azi1 azimax=azi2 if(nazi.gt.1) then dazi=1.*(azi2-azi1)/(nazi-1.) do 5 i=1,nazi azitable(i)=azi1+(i-1)*dazi if(azitable(i).lt.azimin) azimin=azitable(i) if(azitable(i).gt.azimax) azimax=azitable(i) 5 continue endif if(nazi.eq.1) then dazi=0. azitable(1)=azi1 endif if(nazi.lt.0) then do 6 i=1,-nazi read(*,*) azitable(i) 6 continue nazi=-nazi endif c azitable = ring longitude from sub-observer point c testfii=-azitable c sunfii=testfii+dlobs do 10 i=1,nazi obstable(i)=-azitable(i) suntable(i)=obstable(i)+dlobs 10 continue read(*,*) nphot,mx,my read(*,*) allfile read(*,*) parfile read(*,*) outfile read(*,*) epsfile read(*,*) wwwfile read(*,*) method read(*,*) ppar1,ppar2,ppar3 read(*,*) icheck read(*,*) idum read(*,*) sredu,zredu read(*,*) sizemin read(*,*) s20,f20 if(method.eq.1 ) smet='LAMBERT-element ' if(method.eq.11) smet='LAMBERT-sphere-center ' if(method.eq.12) smet='LAMBERT-sphere-point ' if(method.eq.2 ) smet='ISOTROPIC-sphere' if(method.eq.3 ) smet='HG-sphere-center' if(method.eq.4 ) smet='HG-sphere-point' if(method.eq.13) smet='HG2-sphere-center' if(method.eq.14) smet='HG2-sphere-point' if(method.eq.5 ) smet='EULER-sphere-center' if(method.eq.6 ) smet='EULER-sphere-point' if(method.eq.7 ) smet='CALLISTO-sphere-center' if(method.eq.8 ) smet='CALLISTO-sphere-point' if(method.eq.9 ) smet='POWER-sphere-center' if(method.eq.10) smet='POWER-sphere-point' nphot_ori=nphot c----------------------------------------- c Saturn-shine: c iglobe=0 -> SUN c iglobe=1 -> SATURN c iglobe=2 -> SATURN-SOUTH c iglobe=3 -> SATURN-NORTH c abox = calculation distance c rsaturn = Saturn radius c nlat_sat = number of latitude bins c nfii_sat = number of longitude bins c model_sat = saturn model c = 1 -> Lambert c = 2 -> RED (DONES ET ALL) c = 3 -> BLUE (DONES ET ALL) rsaturn=60.4 read(*,*) iglobe read(*,*) abox read(*,*) nlat_sat,nfii_sat read(*,*) model_sat if(model_sat.eq.1) smodel='Lambert-planet' if(model_sat.eq.2) smodel='RED 0.5*zones+0.5*belts' if(model_sat.eq.3) smodel='BLUE 0.5*zones+0.5*belts' if(model_sat.eq.4) smodel='RED+BLUE 0.5*zones+0.5*belts' if(model_sat.eq.5) smodel='ZONES 0.75*blue+0.25*red' if(model_sat.eq.6) smodel='BELTS 0.75*blue+0.25*red' c-------------------------------------------------- c correct obsfii,obslat for geometric distortion? c d_per_a = distance of Voyager from Saturn divided by c distance of ring from Saturn c = 0 -> assume infinity, do not correct read(*,*) d_per_a call tim(cpu0) c*********************************************************************** c OUTPUT TO LOG_FILE c*********************************************************************** write(*,*) '==========================================' write(*,*) 'mcnewd_asyobs.for (04.01.2002)' call date_and_time (date, time, zone, values) write(*,999) values(1),values(2),values(3),values(5),values(6) 999 format(1x,i4,'-',i2,'-',i2,2x,i2,':',i2,':',i2) write(*,*) '==========================================' xlonsun=dlobs xlonobs=0. write(*,*) 'RUNID= ',runide(1:len) write(*,*) 'POSDAT= ',ipos1,ipos2,npos write(*,*) 'AZIRANGE= ',azimin,azimax,nazi write(*,*) 'BSUN, LONSUN=',bsun,xlonsun write(*,*) 'BOBS, LONOBS=',bobs,xlonobs write(*,*) 'NPHOT,MX,MY= ',nphot,mx,my write(*,*) 'METHOD = ',method,' ',smet write(*,*) 'PPAR 1,2,3 = ',ppar1,ppar2,ppar3 write(*,*) 'SEED = ',idum write(*,*) 'allfile,parfile,outfile=' write(*,*) '-> ',allfile write(*,*) '-> ',parfile write(*,*) '-> ',outfile write(*,*) '-> ',epsfile write(*,*) '-> ',wwwfile write(*,*) 'sredu,zredu',sredu,zredu write(*,*) 'sizemin ',sizemin write(*,*) 's20,f20 ',s20,f20 if(iglobe.eq.0) write(*,*) 'ILLUMINATION BY SUN' if(iglobe.eq.1) write(*,*) 'ILLUMINATION BY SATURN' if(iglobe.eq.2) write(*,*) 'ILLUMINATION BY SATURN-SOUTH' if(iglobe.eq.3) write(*,*) 'ILLUMINATION BY SATURN-NORTH' if(iglobe.ge.1) then write(*,*) 'MODEL_SAT=',model_sat,smodel write(*,*) 'DISTANCE =',abox,abox/rsaturn write(*,*) 'NLAT,NFII=',nlat_sat,nfii_sat endif write(*,*) 'D_VGR/ABOX (0-> infinite)',d_per_a write(*,*) '==========================================' c*********************************************************************** c OUTPUT TO ALLFILE c*********************************************************************** close(40) open(unit=40,file=allfile,status='unknown') versio=5 write(40,*) versio write(40,*) runide write(40,*) ipos1,ipos2 write(40,*) bsun,bobs,dlobs write(40,*) azi1,azi2,nazi c versio=0 c write(40,*) nphot c versio=1 write(40,*) nphot,mx,my write(40,*) (azitable(i),i=1,nazi) write(40,*) (obstable(i),i=1,nazi) write(40,*) (suntable(i),i=1,nazi) c versio=2 write(40,*) method write(40,*) sredu,zredu write(40,*) sizemin,s20,f20 c versio=3 write(40,*) ppar1,ppar2,ppar3 c versio=4 write(40,*) iglobe write(40,*) abox,rsaturn write(40,*) nlat_sat,nfii_sat write(40,*) model_sat c versio=5 write(40,*) d_per_a c***************************************************************************** c LOOP OVER POSITION_FILES do 100 ipos=ipos1,ipos2 write(40,*) ipos if(ipos.eq.0) posfile=runide(1:len)//'.posdat00' if(ipos.eq.1) posfile=runide(1:len)//'.posdat01' if(ipos.eq.2) posfile=runide(1:len)//'.posdat02' if(ipos.eq.3) posfile=runide(1:len)//'.posdat03' if(ipos.eq.4) posfile=runide(1:len)//'.posdat04' if(ipos.eq.5) posfile=runide(1:len)//'.posdat05' if(ipos.eq.6) posfile=runide(1:len)//'.posdat06' if(ipos.eq.7) posfile=runide(1:len)//'.posdat07' if(ipos.eq.8) posfile=runide(1:len)//'.posdat08' if(ipos.eq.9) posfile=runide(1:len)//'.posdat09' if(ipos.eq.10) posfile=runide(1:len)//'.posdat10' if(ipos.eq.11) posfile=runide(1:len)//'.posdat11' if(ipos.eq.12) posfile=runide(1:len)//'.posdat12' if(ipos.eq.13) posfile=runide(1:len)//'.posdat13' if(ipos.eq.14) posfile=runide(1:len)//'.posdat14' if(ipos.eq.15) posfile=runide(1:len)//'.posdat15' if(ipos.eq.16) posfile=runide(1:len)//'.posdat16' if(ipos.eq.17) posfile=runide(1:len)//'.posdat17' if(ipos.eq.18) posfile=runide(1:len)//'.posdat18' if(ipos.eq.19) posfile=runide(1:len)//'.posdat19' if(ipos.eq.20) posfile=runide(1:len)//'.posdat20' if(ipos.eq.21) posfile=runide(1:len)//'.posdat21' if(ipos.eq.22) posfile=runide(1:len)//'.posdat22' if(ipos.eq.23) posfile=runide(1:len)//'.posdat23' if(ipos.eq.24) posfile=runide(1:len)//'.posdat24' if(ipos.eq.25) posfile=runide(1:len)//'.posdat25' if(ipos.eq.26) posfile=runide(1:len)//'.posdat26' if(ipos.eq.27) posfile=runide(1:len)//'.posdat27' if(ipos.eq.28) posfile=runide(1:len)//'.posdat28' if(ipos.eq.29) posfile=runide(1:len)//'.posdat29' if(ipos.eq.30) posfile=runide(1:len)//'.posdat30' if(ipos.eq.31) posfile=runide(1:len)//'.posdat31' if(ipos.eq.32) posfile=runide(1:len)//'.posdat32' if(ipos.eq.33) posfile=runide(1:len)//'.posdat33' if(ipos.eq.34) posfile=runide(1:len)//'.posdat34' if(ipos.eq.35) posfile=runide(1:len)//'.posdat35' if(ipos.eq.36) posfile=runide(1:len)//'.posdat36' if(ipos.eq.37) posfile=runide(1:len)//'.posdat37' if(ipos.eq.38) posfile=runide(1:len)//'.posdat38' if(ipos.eq.39) posfile=runide(1:len)//'.posdat39' if(ipos.eq.40) posfile=runide(1:len)//'.posdat40' if(ipos.eq.41) posfile=runide(1:len)//'.posdat41' if(ipos.eq.42) posfile=runide(1:len)//'.posdat42' if(ipos.eq.43) posfile=runide(1:len)//'.posdat43' if(ipos.eq.44) posfile=runide(1:len)//'.posdat44' if(ipos.eq.45) posfile=runide(1:len)//'.posdat45' if(ipos.eq.46) posfile=runide(1:len)//'.posdat46' if(ipos.eq.47) posfile=runide(1:len)//'.posdat47' if(ipos.eq.48) posfile=runide(1:len)//'.posdat48' if(ipos.eq.49) posfile=runide(1:len)//'.posdat49' if(ipos.eq.50) posfile=runide(1:len)//'.posdat50' if(ipos.eq.51) posfile=runide(1:len)//'.posdat51' if(ipos.eq.52) posfile=runide(1:len)//'.posdat52' if(ipos.eq.53) posfile=runide(1:len)//'.posdat53' if(ipos.eq.54) posfile=runide(1:len)//'.posdat54' if(ipos.eq.55) posfile=runide(1:len)//'.posdat55' if(ipos.eq.56) posfile=runide(1:len)//'.posdat56' if(ipos.eq.57) posfile=runide(1:len)//'.posdat57' if(ipos.eq.58) posfile=runide(1:len)//'.posdat58' if(ipos.eq.59) posfile=runide(1:len)//'.posdat59' if(ipos.eq.60) posfile=runide(1:len)//'.posdat60' if(ipos.eq.61) posfile=runide(1:len)//'.posdat61' if(ipos.eq.62) posfile=runide(1:len)//'.posdat62' if(ipos.eq.63) posfile=runide(1:len)//'.posdat63' if(ipos.eq.64) posfile=runide(1:len)//'.posdat64' if(ipos.eq.65) posfile=runide(1:len)//'.posdat65' if(ipos.eq.66) posfile=runide(1:len)//'.posdat66' if(ipos.eq.67) posfile=runide(1:len)//'.posdat67' if(ipos.eq.68) posfile=runide(1:len)//'.posdat68' if(ipos.eq.69) posfile=runide(1:len)//'.posdat69' if(ipos.eq.70) posfile=runide(1:len)//'.posdat70' if(ipos.eq.71) posfile=runide(1:len)//'.posdat71' if(ipos.eq.72) posfile=runide(1:len)//'.posdat72' if(ipos.eq.73) posfile=runide(1:len)//'.posdat73' if(ipos.eq.74) posfile=runide(1:len)//'.posdat74' if(ipos.eq.75) posfile=runide(1:len)//'.posdat75' if(ipos.eq.76) posfile=runide(1:len)//'.posdat76' if(ipos.eq.77) posfile=runide(1:len)//'.posdat77' if(ipos.eq.78) posfile=runide(1:len)//'.posdat78' if(ipos.eq.79) posfile=runide(1:len)//'.posdat79' if(ipos.eq.80) posfile=runide(1:len)//'.posdat80' if(ipos.eq.81) posfile=runide(1:len)//'.posdat81' if(ipos.eq.82) posfile=runide(1:len)//'.posdat82' if(ipos.eq.83) posfile=runide(1:len)//'.posdat83' if(ipos.eq.84) posfile=runide(1:len)//'.posdat84' if(ipos.eq.85) posfile=runide(1:len)//'.posdat85' if(ipos.eq.86) posfile=runide(1:len)//'.posdat86' if(ipos.eq.87) posfile=runide(1:len)//'.posdat87' if(ipos.eq.88) posfile=runide(1:len)//'.posdat88' if(ipos.eq.89) posfile=runide(1:len)//'.posdat89' if(ipos.eq.90) posfile=runide(1:len)//'.posdat90' if(ipos.eq.91) posfile=runide(1:len)//'.posdat91' if(ipos.eq.92) posfile=runide(1:len)//'.posdat92' if(ipos.eq.93) posfile=runide(1:len)//'.posdat93' if(ipos.eq.94) posfile=runide(1:len)//'.posdat94' if(ipos.eq.95) posfile=runide(1:len)//'.posdat95' if(ipos.eq.96) posfile=runide(1:len)//'.posdat96' if(ipos.eq.97) posfile=runide(1:len)//'.posdat97' if(ipos.eq.98) posfile=runide(1:len)//'.posdat98' if(ipos.eq.99) posfile=runide(1:len)//'.posdat99' write(*,*) 'posfile=',posfile c******************************************************************************** c LOOP OVER RING_LONGITUDS do 110 iazi=1,nazi ntest=1 testfii(1)=obstable(iazi) sunfii=suntable(iazi) testlat(1)=bobs sunlat=bsun c Geometric correction? if(d_per_a.ne.0) then xlatvgr=bobs xlonvgr=testfii(1) vx=cos(ar*xlonvgr)*cos(xlatvgr*ar) vy=sin(ar*xlonvgr)*cos(xlatvgr*ar) vz= sin(xlatvgr*ar) dx=vx*d_per_a-1. dy=vy*d_per_a dz=vz*d_per_a dd=sqrt(dx**2+dy**2+dz**2) testfii(1)=atan2(dy/dd,dx/dd)*ra testlat(1)=asin(dz/dd)*ra endif if(versio.ge.5) write(40,*) testfii(1),testlat(1) c Saturn-shine? nphot_tod=nphot_ori flux_sat=0. if(iglobe.ge.1) then call saturn_illumination + (sunlat,sunfii,nphot_ori,nphot_tod,flux_sat) nphot=nphot_tod endif c versio=4 write(40,*) flux_sat,nphot c MC-calculation call mcnew_sub(nphot,posfile,parfile,outfile, + sunfii,sunlat, + ntest,testfii,testlat,testw,testd,testn,testm, + idum,icheck,mx,my, + npart,npartori,npot_mean,nbox_mean,iscat_mean, + method,epsfile,wwwfile,ifor,sredu,tau,ndirect, + zredu,ppar1,ppar2,ppar3) call tim(cpu2) write(*,*) 'pos,azi,w,cpu',ipos,azitable(iazi),testw(1), + cpu2-cpu0 if(iglobe.ge.1) write(*,*) 'flux_sat,nphot',flux_sat,nphot c store results nphot_apu=1 if(nphot_apu.gt.0) nphot_apu=nphot table_fdirect(iazi)=1.*ndirect/nphot_apu tablew(iazi)=testw(1) tabled(iazi)=testd(1) tablen(iazi)=testn(1) do 111 jj=1,10 tablem(iazi,jj)=testm(1,jj) 111 continue c------------------------------------------------------------------ c end loop over ring-longitudes 110 continue write(40,*) (tablew(i),i=1,nazi) write(40,*) (tabled(i),i=1,nazi) write(40,*) (tablen(i),i=1,nazi) write(40,*) ((tablem(i,j),i=1,nazi),j=1,10) c versio=2 write(40,*) (table_fdirect(i),i=1,nazi) c------------------------------------------------------------------ c end loop over position-files 100 continue write(40,*) cpu2-cpu0 write(40,*) npart,npartori,npot_mean,nbox_mean,iscat_mean write(40,*) tau write(*,*) 'npart,npot,nbox,iscat', + npart,npot_mean,nbox_mean,iscat_mean write(*,*) '==========================================' call date_and_time (date, time, zone, values) write(*,999) values(1),values(2),values(3),values(5),values(6) write(*,*) '==========================================' close(40) end